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ABSTRACT 


The  behavior  of  the  electric  field  together  with  the  electron  and  ion  densities  in  the 
vicinity  of  a  non-emitting,  plane  anode  is  investigated.  The  selected  approach  involves  non- 
linear analysis  techniques  on  the  continuum  equations  for  steady-state,  isothermal  conditions 
where  both  ionization  and  two-body  recombination  are  included.  Ions,  created  through 
electron  bombardment  of  neutral  atoms,  are  repelled  toward  two  stagnation  regions:  within 
or  near  the  sheath  boundary  and  near  the  plasma  interface.  These  equilibria  form  as  a  result 
of  the  chemistry  present:  recombination  establishes  the  latter  while  ionization  stipulates  the 
former.  As  presented,  the  sheath  is  fundamentally  unstable  -  ions  are  driven  toward  the 
negative  electrode.  Using  nitrogen  data  for  a  numeric  example,  the  following  observations 
are  made:  a  sufficiently  strong  applied  electric  field  pushes  the  ion  density  toward  that  of  the 
electrons  through  a  well  -  a  constrictive  phenomenon.  Both  a  transition  region,  dominated 
by  density  gradients,  and  a  diffusion-driven  zone  are  found  to  move  the  system  toward  the 
plasma  interface.  The  characteristics  of  this  process  are  influenced  by  the  applied  electric 
field,  but  the  instability  of  the  chemistry-induced  stagnation  regions  precludes  numeric 
convergence.  Insufficient  dissipation  may  prevent  the  stability  of  the  anode  fall  model  as 
presented.  Suggested  improvements  to  the  model  descriptions  include  considering  the  effects 
of  temperature  gradients,  magnetic  fields,  three-body  recombination,  diffusion  written  in 
terms  of  the  electric  field,  multi-dimensionality  and/or  time-dependencies. 
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I.  INTRODUCTION 

Whenever  an  object  starts  moving  in  a  given  direction,  a  second  object  must  be 
present  which  accepts  the  reaction  of  the  force  that  accelerates  the  first,  so  says  Newton's 
third  law  of  motion.  When  a  spacecraft  moves  through  empty  space,  the  accelerating  forces 
act  between  the  vehicle  and  the  mass  of  the  propellant  expelled.  In  1895,  the  distinguished 
Russian  rocket  pioneer  Konstantin  Eduarovitch  Tsiolkovskii  derived  the  famous  "Tsiolkovskii 
equation,"  the  basis  of  all  theoretical  work  on  rocket  propulsion.  Relating  the  burnout 
velocity  of  a  rocket  with  its  exhaust  velocity  and  its  mass  ratio,  this  "rocket  equation"  shows 
how  important  it  is  for  the  exhaust  velocity  of  a  rocket  motor  to  be  as  high  as  possible.  [Ref. 
1:  pp.  2-3] 

Since  mankind's  advent  into  space,  several  types  of  space  flight  propulsion  systems 
have  been  developed  to  include  chemical,  nuclear,  electric  and  solar  propulsion  devices  with 
the  majority  of  spacecraft  to  date  employing  chemical  rockets.  While  the  science  and 
technology  of  chemical  rockets  developed  at  a  rapid  pace,  the  idea  of  using  the  benefits  of 
Maxwell's  relations  to  eject  charged  particles  from  a  thrust  chamber  progressed  slowly  with 
long  intervals  between  periods  of  creative  thought.  In  1945,  Herbert  Radd  recommended  the 
use  of  fast  electrically  accelerated  particles  to  reduce  fuel  mass  and  in  1946,  Jakob  Ackeret 
analyzed  the  performance  capability  of  a  spacecraft  by  arguing  that  the  total  kinetic  energy 
carried  away  by  the  exhaust  beam  corresponded  to  a  reduction  of  the  mass  of  the  vehicle 
according  to  E=mc2.  Hence,  assuming  that  a  certain  fraction  of  the  vehicle  mass  is  converted 
into  kinetic  beam  energy,  the  total  propellant  mass  can  be  calculated  which  in  turn  leads  to 
a  maximum  terminal  velocity  of  the  spacecraft,  a  result  which  is  independent  of  particle  size, 
time  of  propulsion,  conversion  method  and  magnitude  of  the  vehicle.  In  essence,  the 
optimum  exhaust  velocity  is  determined  simply  by  maximizing  the  terminal  velocity.  [Ref.  1 : 
pp.  2-3] 

Propulsive  devices  that  use  Maxwell's  Law's  differ  from  chemical  systems:  the  latter 
is  characterized  by  having  the  same  energy  and  propulsive  sources,  while  in  the  former  the 
energy  conversion  method  differs  from  that  of  propulsion.   In  essence,  electric  propulsion 
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systems  allow  greater  carrying  capacity  with  shorter  flight  times.  For  example,  the  Hohmann 
transfer  time  from  a  Low-Earth-Orbit  (LEO)  to  the  orbit  of  Pluto  requires  approximately  45 
years  if  chemical  means  are  considered;  an  electrically  propelled  vehicle,  in  comparison, 
requires  approximately  three  years.  [Ref.  1:  p.  viii].  The  large  difference  in  mass  payload  ratio 
(final  mass/initial  launch  mass)  obtained  is  highlighted  by  an  analysis  for  a  Mars  mission:  a 
chemical  system  using  a  high  impulse  Hohmann  trajectory  from  LEO  delivers  approximately 
10%  to  18%  of  the  launch  mass  to  the  Martian  surface  [Ref.  2:  pp.  152],  while  an  electric 
system  using  a  low  impulse  spiral  trajectory  could  deliver  20%  to  60%  of  the  launch  mass, 
depending  on  the  desired  transit  time.  With  effective  exhaust  velocities  as  high  as  10,000  m/s, 
electric  propulsion  thus  offers  the  performance  envelope  needed  for  manned  interplanetary 
missions.  [Ref.  2:  pp.  37] 

Electric  propulsion  systems  are  divided  into  three  categories:  1)  electrostatic,  where 
acceleration  is  achieved  through  the  interaction  of  coulombic  fields  with  charged  propellant 
particles,  2)  electrothermal,  where  propellant  heated  electrically  is  expanded 
thermodynamically  through  a  nozzle,  and  3)  electromagnetic,  where  acceleration  is  achieved 
through  the  interaction  of  the  Lorentz  force  with  highly  ionized  gases  (plasma).  Of  the  types 
of  electric  propulsive  thrusters  available,  electrothermal  devices  include  the  resistojet  and 
arcjet.  It  should  be  noted  that  the  current  necessary  to  form  the  plasma  in  these  devices  also 
allows  the  magnetic  field  to  accelerate  the  gas.  These  systems  (and  for  that  matter,  any 
device  using  a  thermal  nozzle)  are,  however,  limited  by  the  amount  of  heat  that  can  be  added 
to  the  flow  and  are  thus  ultimately  dependent  on  the  materials  used.  As  electromagnetic 
(EM)  thrusters  rely  on  Lorentz  interactions  to  accelerate  the  flow,  these  devices  are  not 
limited  in  the  sense  of  electrothermal  systems  and  therefore  offer  different  operational 
characteristics.  The  various  types  of  EM  thrusters  available  include  the  ion  engine  and  plasma 
thrusters  (wide  variety  to  include  electrodynamic,  magnetoplasma-dynamic  or  MPD,  pulsed- 
plasma,  Hall  accelerator,  Lorentz  force  accelerator,  pulsed  coaxial  thruster,  plasma  arcjet, 
pulsed  ablative  thrusters,  to  name  just  a  few).  [Ref.  2:  pp.  565-596]  The  MPD  thruster  is 
relevant  to  this  work. 

MPD  systems  are  a  relatively  recent  addition  to  the  electric  propulsion  family, 


originally  discovered  when  peculiar  behavior  of  an  electrothermal  arcjet  thruster  was 
investigated.  Researchers  found  that  without  significant  mass  flow-rate  applied,  the  arcjet 
continued  to  operate  in  a  mode  clearly  unique:  pressure  within  the  thrust  chamber  was  too 
low  for  the  device  to  be  choked,  and  so  it  was  theorized  that  an  induced  magnetic  field  from 
the  high  arc  current  accelerated  the  plasma  flow  -  the  magnetoplasmadynamic  thruster  was 
"born."  [Ref.  3] 

The  area  of  research  to  which  MPD  thrusters  are  applied  is  termed 
magnetohydrodynamics  (MHD),  the  field  of  physics  that  investigates  how  moving  conductive 
(electrolytic)  fluids  behave  while  subjected  to  applied  electric  or  magnetic  fields.  The  fluids 
involved  must  be  conductive,  or  ionized  sufficiently  to  pass  a  current,  and  in  the  case  of  an 
MPD,  the  fluid  must  reach  a  plasma  state.  MHD  devices  behave  analogously  to  motors  or 
generators  and  can  either  add  or  remove  energy  from  the  fluid  in  exchange  for  electrical 
energy;  systems  developed  include  the  MHD  tunnel  drive  for  submarines  and  the  MHD 
turbine.  Magnetoplasmadynamic  (MPD)  thrusters  are  currently  being  considered  for 
application  to  primary  propulsion  system  on  robotic  and  piloted  interplanetary  missions;  the 
high  specific  impulses  offered  by  MPD  devices  can  reduce  launch  mass  requirements  by  over 
a  factor  of  two  compared  to  chemical  system  if  the  power  to  thrust  conversion  efficiency  is 
over  50%.  [Ref.  4] 

MPD  thrusters  are  classified  by  considering  two  major  features:  1)  the  power  supply 
which  defines  the  continuous,  quasi-steady  and  pulsed  MPD  systems,  and  2)  the  magnetic 
field  which  characterizes  the  self-field  and  applied-field  MPD  thrusters  [Ref.  3].  A 
continuous,  or  steady-state,  thruster  is  one  in  which  the  propellant  flow  and  current 
throughout  the  device  are  continuous,  thereby  allowing  the  system  to  reach  a  high 
temperature  equilibrium  state.  The  quasi-steady  thruster  does  not  reach  equilibria  since  it 
uses  a  continuous  gas  flow  with  rapidly  pulsed  current,  through  which  high  amounts  of  energy 
are  provided  in  short  bursts.  The  pulsed  thruster  is  similar  to  the  quasi-steady  MPD  system 
except  it  uses  less  frequent  pulses;  equilibrium  flow  is  also  not  reached  and  consequently  low 
temperature  materials  may  be  used  for  constructing  such  a  device.  The  applied-  and  self-field 
thrusters  define  the  way  in  which  the  primary  accelerating  field  (magnetic)  is  affected.  The 


applied-field  system  employs  a  large  current  in  an  external  coil,  or  large  permanent  magnets, 
to  accelerate  the  plasma  whereas  in  the  self-field  device  the  current  not  only  ionizes  the  gas 
but  also  creates  the  field  that  accelerates  the  resulting  plasma. 

The  magnetoplasmadynamic  thruster  was  originally  discovered  in  the  1960's,  and  yet 
much  work  remains  to  understand  its  fundamental  operation.  For  example,  the  steady  state, 
self-field  MPD  thruster  operated  by  the  University  of  Southampton,  U.K.,  was  originally 
developed  to  simulate  the  atmosphere  effects  on  structures  in  LEO,  an  environment 
dominated  by  highly  reactive  atomic  oxygen.  The  MPD  is  suited  to  simulate  such  conditions 
as  it  offers  a  unique  combination  of  high  thrust  and  high  exit  velocity.  The  Southampton 
MPD  presently  operates  at  power  levels  in  the  range  6-17  kW  but  surfers  from  electrode 
erosion,  instabilities,  a  contaminated  beam  and  short  run-times.  Work  on  this  thruster  has 
identified  electrode  erosion  as  a  major  problem,  making  the  device  unsuitable  for  its  intended 
function.  Specifically,  the  mode  of  erosion  differs  from  that  seen  at  other  institutions  and  the 
levels  of  erosion  are  higher  and  exhibit  a  rate  increase  not  seen  with  other  devices.  [Ref.  5] 
Electrode  erosion,  current  spotting,  frozen  flow  loses  and  electrode  power  deposition  are  just 
several  factors  that  limit  the  performance  of  MHD  devices  [Ref.  6,  7].  Specifically,  anode 
power  deposition  is  the  single  largest  power  loss  mechanism  in  MPD  thrusters  operating  at 
sub-megawatt  power  levels.  [Ref.  8,  9,  10,  11,  12] 

To  gain  further  insight  to  MPD  phenomena,  computer  codes  that  accurately  describe 
observed  data  from  steady-state  MPD  thrusters  have  been  developed  [Ref.  13,  14,  15]. 
However,  these  codes  do  not  adequately  describe  observed  data  from  quasi-steady  thruster 
experiments.  It  has  been  suggested  that  the  lack  of  proper  electrode  modeling  (i.e.,  sheaths 
and  fall  potentials)  in  these  codes  may  explain  this  discrepancy  [Ref.  16];  perhaps  the  difficult 
set  of  coupled,  nonlinear  partial  differential  equations  involved  has  limited  analytical  work  to 
model  the  anode  sheath  and  ambipolar  regions.  For  example,  Hugel  [Ref.  17]  and 
Subramanian  [Ref.  18]  address  the  influence  of  the  sheath  region,  but  do  not  model  the 
electric  field,  temperature  or  sheath  fall  voltage.  Hence,  the  focus  of  this  work:  analysis  of 
the  anode  voltage  drop  region. 


H.  LITERATURE  REVIEW 

The  principal  voltage  loss  mechanism  near  an  electrode  can  be  divided  into  two  main 
classes:  ohmic  and  sheath  losses.  Ohmic  resistivity  occurs  because  of  the  finite  conductivity 
of  a  real  plasma,  particularly  thermal  boundary  layers,  degree  and  kinetics  of  ionization  and 
Joule  heating.  Sheath  drops  are  caused  by  Debye  shielding  which  forms  a  space  charge  field 
adjacent  to  the  electrode  [Ref.  19:  pp.  200-204].  Furthermore,  material  problems  restrict  the 
temperatures  at  which  the  electrodes  can  operate.  In  many  cases  cooling  of  the  electrodes 
is  required  since  the  plasma,  in  order  to  maintain  a  high  ionization,  must  be  hotter  than  the 
working  temperatures  of  most  materials.  This  temperature  difference  between  the  electrodes 
and  the  plasma  further  aggravates  the  voltage  losses  because  of  the  presence  of  the  thermal 
boundary  layers.  As  has  been  shown,  voltage  losses  can  exceed  50%  of  the  total  power 
output,  with  sheath  drops  accounting  for  a  large  fraction  of  this  loss.    [Ref.  6,  7] 

The  tool  most  commonly  used  for  measurement  of  local  plasma  properties  is  the 
electrostatic  probe,  a  small  electrode  which  may  be  biased  positively  or  negatively  with 
respect  to  the  plasma  in  which  it  is  immersed.  Such  work  is  relevant  to  non-emitting  anodes 
which,  though  heavily  biased,  also  probe  and  disturb  the  plasma  locally.  As  a  result,  most 
analytical  work  on  sheath  phenomena  is  embodied  in  probe  theory  investigations  where  the 
effects  of  an  electrode  (probe)  locally  disturb  a  quiescent  plasma.  [Ref.  20,21:  pp.  1-7,  pp. 
8-9] 

In  1 924  Langmuir  pioneered  the  use  of  electric  probes,  giving  the  first  sound 
mathematical  basis  to  this  diagnostic  technique.  The  analysis  of  electrostatic  probes  in 
collision-dominated  (higher  pressure)  plasmas  began  with  Davydov  and  Zmanovskaja  [Ref. 
22]  who  assumed  that  a  quasi-neutral  diffusion  controlled  region  extended  to  within  one  mean 
free  path  of  the  probe  where  it  was  matched  to  the  edge  of  a  free-fall  sheath;  a  transition 
region  between  the  two  regimes  was  not  attempted,  however.  Ecker,  et  al.  [Ref.  23] 
attempted  matching  the  solution  in  the  collision-dominated  region  directly  to  that  in  the  free- 
fall  sheath  using  variational  principles,  work  which  was  criticized  by  Cohen  [Ref.  24]  based 
on  mathematical  grounds:  something  was  forced  to  fit  via  some  artificially  imposed  boundary; 


matching  between  separate  regions  was  accomplished  through  variational  principles  and 
consequently  was  not  asymptotic  in  nature  [Ref.  24]. 

Existing  sheath  solutions  are  limited  to  special  case  simplifications.  Several  solutions 
are  available  for  spherical  probes  by  Kiel  [Ref.  25],  whose  analysis  neglected  diffusion,  Barad 
and  Cohen  [Ref.  26],  Su  and  Lam  [Ref.  27],  as  well  as  by  Cohen  [Ref  24,  28].  The  latter 
group  carried  out  the  first  completely  systematic  analyses  of  spherical  probes  in  collision- 
dominated  plasmas.  These  researchers  assumed  the  continuum-equations  valid  everywhere 
in  the  plasma  right  up  to  the  probe  surface;  their  results  showed  both  the  sheath  and  the  quasi- 
neutral  regions  appeared  naturally  from  the  diffusion  equations.  In  these  efforts,  two  major 
assumptions  prevail:  the  plasma  was  assumed  so  lightly  ionized  that  only  charge-neutral 
interactions  were  accounted,  and  the  continuum  equations  were  postulated  with  an  isothermal 
formulation.  For  cool,  nonequilibrium  plasmas,  however,  charge-charge  collisions  must  be 
considered  [Ref.  29:  p.  165].  To  that  end,  McKee  and  Mitchner  [Ref.  30]  investigated  a 
collisionless  sheath  with  constant  transport  properties  that  also  included  ionization  as  well  as 
recombination  effects  in  the  ambipolar  region.  Exploring  spherical  probe  theory  further, 
Cohen  and  Schweitzer  [Ref.  31]  then  extended  Cohen's  original  asymptotic  analysis  [Ref.  24] 
for  the  entire  region,  an  analysis  which  allowed  for  weak  ionization  and  recombination  effects 
while  also  assuming  constant  transport  coefficients.  Su  and  Sonin  [Ref.  32]  next  formulated 
a  spherical  probe  theory  by  considering  the  effects  of  electron-ion  collisions  by  approximating 
diffusion  coefficients  as  a  function  of  the  neutral  particle  density,  a  density  assumed  as 
constant.  Barad  and  Cohen  [Ref.  26]  then  presented  a  continuum  theory  (for  spherical 
probes)  in  which  non-constant  transport  coefficients  were  used,  coefficients  based  on  charge- 
charge  collisional  assumptions.  This  analysis,  however,  did  not  incorporate  ionization  or 
recombination  effects.  In  work  similar  to  Lam's  original  theory  for  flow  of  weakly  ionized 
gases  [Ref.  28],  Stahl  and  Su  [Ref.  33]  use  the  same  approach  of  separating  the  sheath, 
ambipolar  region  and  free  stream;  their  methods  prove  the  existence  of  a  sheath  on  a  flat 
probe,  yet  plasma  chemistry  is  not  considered  (ionization  and  recombination).  It  should  be 
noted  that  both  Su  and  Lam  [Ref.  27],  for  spherical  probes,  and  Stahl  and  Su  [Ref  33],  for 
flat  probes,  consider  transition  regions  to  match  the  inner  sheath  to  the  outer  quasi-neutral 


or  ambipolar  zones  -  albeit  both  research  groups  neglected  plasma  reactivity  in  their 
formulations.  In  related  work,  Godyak  and  Sternberg  [Ref.  34]  applied  the  ionization 
potential  as  the  boundary  condition  for  both  the  plasma  and  the  sheath  in  an  effort  to 
understand  the  electrodynamic  properties  of  the  plasma-sheath  boundary  for  the  cathode;  their 
analysis,  however,  assumed  constant  ionization  for  the  plasma  reactivity  with  recombination 
neglected. 

It  may  be  inferred  from  these  efforts  that  anode-voltage-drop-region  analyses  remain 
insoluble  due  not  only  to  the  nonlinear  nature  of  the  equations,  but  also  to  the  transitional 
nature  of  particle  motion  in  most  sheaths.  In  an  extensive  review  of  probe  theory,  Chung, 
Talbot  and  Touryan  [Ref.  20]  state  that  a  general  solution  for  determining  charge  density  and 
species  temperature  for  probes  small  relative  to  boundary  layer  thicknesses  is  not  available. 
Additionally,  Nasser  [Ref.  35]  discusses  an  elementary  theoretical  approach  to  the  glow 
discharge  problem  by  introducing  a  set  of  one-dimensional  governing  equations  (continuum)  - 
the  electron  and  ion  current  equations,  the  net  production  of  electrons  and  ions,  and  the 
Poisson  equation.  According  to  Nasser,  these  equations  should  have  a  solution,  but  most 
attempts  to  achieve  such  have  failed,  with  the  difficulty  lying  in  the  boundary  conditions  [Ref. 
35:  p.  412].  Voltage  losses  at  electrode  boundaries  and  surface  erosion,  together  with  sheath 
effects  which  influence  ionization,  play  an  important  role  in  plasma  devices  and  must  be 
controlled  in  designs  of  practical  interest.  Particularly,  thermal  arcjets  and  MPD  accelerators 
deposit  between  15%  and  80%  of  the  input  power  into  the  anode.  This  presents  not  only  a 
severe  performance  penalty,  but  also  introduces  thermal  design  problems  since  the  heat  thus 
generated  must  be  radiated  away  from  the  thruster  surfaces  [Ref.  7,  9]. 

Completely  satisfactory  solutions  to  the  anode  voltage  drop  phenomena  do  not  exist 
for  several  plasma  conditions,  one  such  case  involves  steady,  collisional,  low-temperature 
isothermal  discharges  [Ref.  36].  For  instance,  the  effect  of  several  parameters  on  the  region 
require  further  study,  including  the  influence  of  temperature,  the  role  of  the  problem 
dimensionality  as  well  as  the  effect  of  plasma  chemistry.  It  remains  unclear,  moreover,  as  to 
the  exact  cause  of  the  anode  spots  which  have  been  observed  for  various  MPD  thruster 
operating  conditions,  a  behavior  which  appears  roughly  as  periodic  oscillations  [Ref.  6,  7,  9]. 


A  one-dimensional,  collisional,  equilibrium  solution  can  satisfactorily  reproduce  the  observed 
electric  field  and  charge  density  distributions  for  the  entire  anode  voltage  drop  region  where 
the  electron  temperature  equals  that  of  the  heavy  species  [Ref.  36].  However,  this  model 
cannot  describe  any  decrease  in  current  density  away  from  the  surface,  or  current  constriction, 
at  the  anode  surface  which  might  be  necessary  in  non-equilibrium.  A  two-dimensional  model, 
developed  by  Biblarz  and  Dolson  [Ref.  6],  represents  these  phenomena,  in  the  absence  of 
species  ionization  and  recombination,  and  predicts  the  voltage  drop  in  the  region:  the  sheath 
accounted  for  a  majority  of  the  anode  voltage  drop  [Refs.  6,  7,  36]. 

Since  the  pioneering  work  of  Tonks  and  Langmuir  (relating  to  low-pressure 
discharges),  the  anode  description  has  remained  as  one  of  the  classic  problems  in  plasma 
physics.  In  this  work,  the  high  density  problem  is  generalized  to  any  electrode  surface  in 
contact  with  a  partially-ionized  plasma  when  the  current  flow  is  sufficiently  small:  the 
problem  of  a  non-emitting  electrode  (probe)  in  contact  with  a  thermally  generated  plasma. 
An  anode  is  explicitly  considered  along  with  a  non-emitting  cathode  and  the  governing  one- 
dimensional  continuum  equations  are  formulated.  The  electrodynamic  behavior  of  plasma 
reactivity  in  the  electrode  regions  is  then  discussed  with  considerations  given  to  three  distinct 
zones  of  the  anode  voltage  drop  regime:  inner  or  space-charged  sheath,  transition  and  outer 
or  quasi-neutral  ambipolar  zones. 


m.  GOVERNING  EQUATIONS  AND  NUMERICAL  DATA 

A.         GOVERNING  EQUATIONS 

A  flat  electrostatic  probe  immersed  in  a  high-density,  flowing,  partially-ionized  gas  is 
considered.  For  plasmas  of  moderately  high  density  and  low  degree  of  ionization  (as  applied 
to  discharge  lasers  and  MHD  devices),  the  plasma  region  affected  by  an  electrode,  i.e.,  the 
electrode  region,  can  be  considered  of  the  boundary-layer  type  with  the  transport  of  charges 
approximated  as  one-dimensional.  Traditionally,  the  electrode  region  is  further  separated  into 
two  distinct  zones,  as  delineated  by  the  species  densities:  a  sheath  and  an  ambipolar  region, 
both  which  may  be  of  the  order  of  tenths  or  even  hundredths  of  a  millimeter  in  extent  [Ref. 
6]  where  the  sheath  extent  is  driven  by  variation  of  the  electric  field  and  the  quasi-neutral 
region  extent  is  influenced  by  the  influences  of  chemistry. 

The  problem  is  formulated  by  considering  the  hydrodynamic  continuum  equations  for 
two  species,  ions  and  electrons.  Induced  and  applied  magnetic  fields  are  neglected  with 
energy  conservation  implicitly  satisfied  since  heat  transfer,  radiation,  Joule  heating  and  surface 
emission  effects  are  at  first  neglected  (for  example,  B=0,  VTe=0,  VTj=0,  J-E=0). 

The  equations  below  are  assumed  valid  throughout  the  entire  domain  of  interest,  to 
include  the  electrode  (probe)  surface.  Table  (1)  lists  the  nomenclature  used: 

d(Mn.v.) 
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Vj i  =  gV  ( n.  v. )  Equation  (3a) 

V j e  =  qV(neve)  Equation  (3b) 

V  (  n i  v . )  =  v  ne  -  a2  n i  ne  Equation  (4a) 

V(neve)  =-vne  +  cx2n.ne  Equation(4b) 
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Equations  (la)  and  (lb)  describe  the  momentum-transfer  between  the  ions  and  electrons  with 
neutrals:  the  left-hand  side  represents  the  unsteady  and  convective  contributions, 
respectively,  while  the  right-hand  side  symbolizes  coulombic,  pressure,  collisional  and 
frictional  forces,  respectively.  Next,  the  electrode  is  considered  as  a  positive,  non-emitting 
surface  (anode).  Through  this  process,  a  flow  of  negative  charges  is  induced  toward  the  plate 
and  a  flow  of  positive  charges  is  repelled  away  from  the  plate;  the  anode,  in  completing  the 
arc  current  with  the  cathode,  collects  incoming  electrons  while  driving  all  positive  charges 
toward  the  negative  surface.  This  separation  of  charge  in  the  anode  fall  region  produces  a 
divergence  of  the  electric  field  and  a  variation  of  potential  away  from  the  plate,  as  given  by 
equation  (2)  (the  variation  of  potential  is  readily  deduced  from  Gauss's  equation  (2)).  The 
current  densities  are  denoted  by  equations  (3  a)  and  (3  b)  in  terms  of  the  species  flux  of 
charged  particles  per  unit  area.  Equations  (4a)  and  (4b)  are  then  the  continuity  equations  for 
the  ion  and  electron  flux  as  caused  by  the  ionization  and  electron-ion  (two-body) 
recombination  processes:  positive  ions  are  produced  within  the  anode  fall  region  by  electron 
impact  with  neutrals  and  conversely,  positive  ions  combine  with  negative  electrons  to  form 
neutrals.  Equation  (5)  is  the  conservation  of  total  charge;  the  sum  of  the  species  currents 
remains  constant  for  a  one-dimensional  formulation.  Hence,  from  set  (1-5),  the  anode  model 
presented  in  this  work  is  developed  with  the  following  assumptions: 

1 .  Steady  state;  free-stream  plasma  is  neutral  with  uniform  charge  densities  (Sana 
equilibrium). 

2.  Mean  free  path  of  the  charged  particles  is  small  compared  to  the  sheath  thickness 
which  in  turn  is  much  smaller  than  the  dimensions  of  the  body;  the  sheath  is  treated 
as  a  continuum. 

3.  Chemical  equilibrium  in  the  boundary  layer;  the  two-body  recombination 
coefficient  is  held  constant  over  the  entire  anode  fall  region  as  a  first  approximation 
[Ref  29:  p.  238].  Further,  the  ionization  coefficient  is  treated  as  Townsend's  first 
coefficient. 

4.  There  is  no  applied  magnetic  field  and  induced  magnetic  field  effects  are  neglected 
(low  Magnetic  Reynolds  number)  -  consequently  no  ion  slip. 
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5.  No  continuum  radiation  losses,  no  ion  emission,  no  significant  Joule  heating. 

6.  Isothermal:  VT=VTf=VTo=0  and  Tc=Ti=T0.  The  true  electron  temperature  will 
either  be  close  to  the  anode  surface  temperature  or  the  free  stream  temperature  since 
thermal  equilibrium  is  assumed  in  the  boundary  layer.  Because  of  the  relatively  small 
number  of  collisions  that  the  electrons  suffer  in  the  boundary  layer,  it  is  frequently 
argued  that  the  electron  temperature  is  close  to  the  free  stream  temperature  [Ref.  37]. 

7.  Diffusion  velocities  of  the  ions  and  electrons  due  to  electric  field  interactions  are 
small  in  comparison  with  their  thermal  velocities;  diffusion  coefficients  are  assumed 
constant. 

8.  End  effects,  in  the  case  of  a  planar  surface,  are  neglected. 

9.  Under  the  assumption  of  a  highly  biased  probe,  it  is  possible  to  ignore  convection 
not  only  in  the  very  thin  transition  region,  but  also  in  the  sheath  region,  even  for 
sheaths  with  the  same  order  of  thickness  as  the  boundary  layer  [Ref.  33].  Assuming 
perfectly  elastic  collision  between  the  species  and  neutrals,  essentially  a  constant 
collisional  frequency  is  assumed  sufficiently  large  so  that  the  convective  derivative  in 
(1)  can  be  neglected. 

10.  With  the  assumptions  of  steady  state  and  negligible  convection,  a  spatially 
homogenous  fluid  results. 

For  a  semi-infinite  plane  with  coordinate  outward  from  the  planar  positive  surface  so 
that  V=d/dy,  set  (1-5)  thus  becomes: 
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Here,  equations  (3a)  and  (3b)  have  been  incorporated  with  equations  (4a)  and  (4b).  Using 


kT 

species  mobilities,  \x.ie  =  — 2 —    and  the  diffusion  coefficients,  Die  =  — — 

'        miji,e  '        miji,, 


(the  Einstein 


relation   n_  e  =  — l—  ),  species  momentum  conservation  can  be  re- written  in  terms  of  species 

ft  1  -  _ 


flux,  Tie,  where  Fick's  law  of  diffusion  is  a  special  case  (E=0  so  that  Ui,c=0).  Specifically,  (1) 
and  (2)  become 


r      =n. 


ife=ni,evi,e=±VirenireE-D.fVn.ie       Equation  (6) 


Essentially,  diffusion  is  a  random  process  in  which  a  net  flux  from  dense  regions  to  less  dense 
regions  occurs  simply  because  more  particles  begin  their  motion  in  the  more  dense  area  -  the 
flux  is  proportional  to  the  density  gradient.  [Ref  38:  pp.  135-13&]  But  as  the  species  travel 
toward  the  electrode,  reactivity  occurs  -  ions  and  electrons  recombine  while  electron  collision 
with  neutrals  breaks  neutrality  adding  charges  to  the  electrode  region;  flow  neutrality  is 
broken  and  the  fluid  is  no  longer  a  plasma.  Thus,  under  constant-property  conditions,  the 
species  continuity  equations  together  with  Gauss's  equation  describe  the  flux  of  charges  in 
the  vicinity  of  an  electrode  or  plasma  probe.  [Ref.  33,  39]  With  the  application  of  the  Einstein 
relation,  set  (1-5)  are  re-cast  to  give  the  model  used  in  this  work,  the  geometry  for  which  is 
shown  in  Figure  (1): 
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Figure  1:  Geometry  of  Anode  Fall  Region 
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To  obtain  the  mathematical  model  used  in  this  work,  the  definition  of  the  diffusion  coefficient 
is  used  to  re-write  the  species  elastic  collision  frequencies, /s=[kT0]/[Dsms]  for  isothermal 
conditions,  and  the  current  density  is  defined  as  j^r^y.  where  s=i,  e.  The  resulting  equations 
(1-5)  are  then  solved  for  their  highest  derivative  terms: 
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Equation  (7a) 
Equation  (7b) 
Equation  (8) 

Equation  (9a) 
Equation  (9b) 


In  the  formulation  of  equations  (7-9),  14  variables  result  (three  universal  constants 
included).  Before  analysis  of  this  system,  the  equations  must  be  written  in  terms  of 
dimensionless  variables  and  proper  parameters.  For  this  work,  both  the  Buckingham  Pi 
Theorem  and  Fractional  Analysis  techniques  were  applied  (Appendix  A).  The  Buckingham 
method,  although  a  useful  tool  for  revealing  the  dimensionality  of  the  problem,  does  not  bring 
forth  the  physical  effects  which  influence  system  behavior.  On  the  other  hand,  fractional 
analysis  coupled  with  the  use  of  the  governing  equations,  begins  to  give  physical  meaning  to 
the  dimensionless  parameters  obtained. 

Using  undisturbed  plasma  values  for  normalizing  the  variables  gives  three  primary 
length       scales      for      the       system       modeled      by       equations       (7-9)       since 
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.  The  results  are  summarized  in  Table  2  with 
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the  length  scales  listed  in  Table  3: 


Governing  Equations 
[Equations  (7-9)] 
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Table  2:  Non-dimensional  Parameters 
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Table  3:  Length  Scales  via  Fractional  Analysis 
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It  should  be  noted  that  the  length  scales  Lg,  Ln  and  Lj  are  in  agreement  with  reference  50  if 
the  species  conduction  currents  in  the  undisturbed  plasma  are  written  in  terms  of  thermal 
velocities  and  the  above  electric  field  and  number  density  scales  are  combined  to  form  the 
Debye  shielding  length  for  planar  surfaces  ( l^L^X^'2 ). 
B.         NUMERICAL  DATA 

To  gain  further  insight  into  the  physical  models  presented,  numerical  analyses  were 
conducted  with  data  from  reference  36  for  6000°K  nitrogen;  data  for  the  ionization  and  two- 
body  recombination  coefficients  are  found  in  Appendix  B. 


Data  Value 

Temperature  (°K),  T0 

6.0000e+03 

Total  Current  (mA/cm2),  J:  used  for  j^.  j^ 

l.OOOOe+02 

Undisturbed  Plasma  Charge  Density  (1/m3),  n,,:  used  for  n^,  n^ 

1.0000e+19 

Undisturbed  Plasma  Total  Density  (1/m3),  N: 

2.0000e+24 

Undisturbed  Plasma  Electric  Field  (V/m),  E0 

1.2000e+04 

E/N  (Vcm2) 

6.0000e-17 

Pressure  (N/m2),  p:  from  [p=NkT] 

1.6560e+05 

Two-body  recombination  coefficient  (m3/s),  a2:  [Ref.  40] 

1.0000e-13 

Ionization  coefficient  (1/s),  v:  from  [Ref.  41  ]  with  v/N=  10'18  (m3/s)  for  above 
E/N  converted  to  Td  (lTd=10-21Vm2) 

1.0000e+06 

Ion  Diffusion  Coefficient  (m2/s):  D;N=4.46(1018)To1/2  from  [Ref.  42:  pp.  2-28] 

1.7274e-04 

Electron  Diffusion  Coefficient  (m2/s):  D/D^m/M 

3.1665e-01 

Table  4:  Numerical  Data  for  Analysis  -  Nitrogen 
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IV.  PHYSICAL  MODEL 

A.         SPACE-CHARGE  REGION:  SHEATH 

Voltage  losses  at  electrode  boundaries,  surface  erosion  and  sheath  effects  play  an 
important  role  in  plasma  devices  and  must  be  controlled  in  designs  of  practical  interest. 
Particularly,  thermal  arcjets  and  MPD  accelerators  deposit  between  15%  and  80%  of  the 
input  power  into  the  anode.  This  presents  not  only  a  severe  performance  penalty,  but  also 
introduces  thermal  design  problems  since  the  heat  thus  generated  must  be  radiated  away  from 
the  thruster  surfaces.    [Ref.  7,  9]  So  a  natural  question:  what  happens  at  the  anode? 

Near  a  surface  (anode)  in  contact  with  a  plasma  there  is  a  region,  the  sheath,  in  which 
the  electric  field  is  strong  and  in  which  the  electron  number  density  is  much  larger  than  that 
of  the  ions  (see  Figures  (2)  and  (3)). 
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Figure  2:  Axial  Potential  Profile  for  a  High 
Current  Arc  Between  Two  Electrodes  (not 
to  scale)  [after  Ref.  3:  p.  1 14] 
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Figure  3:  Qualitative  Spatial 
Distribution  of  Species  Densities  - 
Anode  Fall  Region  (not  to  scale) 


Figure  (3)  shows  the  qualitative  behavior,  in  the  anode  fall  region,  of  the  spatial 
distribution  of  species  densities;  it  is  expected  that  from  the  wall  outward,  the  densities 
approach  one-another,  eventually  reaching  Sana  equilibrium  at  the  plasma  boundary.   The 
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overall  shape  of  the  curves  (Figure  (3))  stem  purely  from  the  density  flux   —  =  \n  -  a2n 


in  the  absence  of  any  other  types  of  influences.  However,  it  is  usually  the  interaction  of  a 
systems'  many  parts,  not  the  action  of  each  piece  independently,  that  determines  the  final 
system  behavior. 

The  mathematical  model  developed  for  the  sheath,  whether  anode  or  cathode,  depends 
on  the  details  desired:  whether  surface  effects  are  considered,  whether  induced  or  external 
magnetic  fields  are  included,  as  to  how  the  flow  chemistry  is  modeled,  among  others.  For  the 
cathode,  sheath  models  usually  adopt  the  assumptions  that  electrons  are  in  equilibrium  with 
the  electric  field  and  that  the  ionization  in  the  sheath  can  be  neglected  [Ref.  34].  In 
comparison,  the  electrons  in  the  anode  fall  must  supply  by  ionization  enough  ions  to  account 
for  the  ion  current  that  flows  out  of  the  positive  column  toward  the  cathode  [Ref.  43:  pp.  59- 
60];  electrons  emerge  from  the  column,  are  attracted  by  the  anode,  gain  energy  through  the 
voltage  drop  and  ionize  neutrals.  The  extent  of  the  anode  sheath  is  then  determined  by  the 
space  charge  distribution  and  by  the  magnitude  of  the  gas  ionization  potential  [Ref  35:  p. 
420];  hence,  the  associated  length  scale  for  the  sheath  is  governed  by  that  which  results  from 
the  dimensional  analysis  of  the  Gauss  equation  (8). 

Moreover,  after  a  current  flow  is  established  between  the  plates,  the  anode  region  may 
exist  in  a  vapor  that  issues  from  the  electrodes.  In  vacuum  arcs,  Miller  [Ref.  44]  characterizes 
the  anode  region  as  operating  in  one  of  five  distinct  modes,  ranging  from  a  passive,  low 
current  mode  to  a  high  current,  fully  developed  spot  mode.  These  observed  instability 
phenomena  are  similar  to  those  noted  in  MPD  thrusters  -  anode  spot  formation  at  high  current 
clearly  limits  anode  lifetime  [Ref.  6,  7,  9,  17].  What  could  cause  such  "spotting?" 

In  a  comprehensive  review  of  anode-spot  phenomena,  Miller  suggests  [Ref.  44]  that 
a  sudden  increase  in  the  voltage  noise  frequently  indicates  that  an  anode  foot-print  could  be 
forming  or  that  a  transition  into  the  anode-spot  mode  could  be  occurring.  Further,  this  abrupt 
change  in  voltage  has  been  associated  with  a  sudden  change  of  ion  density  in  the  region  near 
the  anode,  a  density  decrease  that  would  leave  the  local  electron  space  charge  uncompensated 
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and  thus  produce  a  low-conductivity  region.  The  sudden  increase  in  voltage  would  then 
reflect  this  local  shortage  of  ions  -  an  ion  depletion  region  or  ion  starvation  region.  Even- 
though  the  most  suitable  model  for  the  anode  region  is  still  actively  debated,  the  general  idea 
of  the  appearance  of  an  ion  deficiency  region  near  the  anode  as  triggering  the  transition  into 
the  anode-spot  mode  has  been  fairly  supported.  [Ref  44]  So,  to  better  understand  the 
processes  involved  within  the  anode  fall,  equations  (7-9)  are  modified  as  follows: 

Because  the  anode-spot  instability  is  likely  related  to  the  species  density  distribution 
within  the  anode  fall  region,  the  role  of  net  charge  production,  the  chemistry,  should  be 
critical  to  the  overall  system  behavior:  from  equations  (7-9),  the  species  densities  are 
influenced  by  both  coulombic  interactions  as  well  as  by  the  currents  and  consequent 
production  phenomena,  the  fine  balance  between  ionization  and  recombination.  Therefore, 
the  current  densities  as  given  by  equations  (9a-9b)  are  combined  with  the  species  densities, 
equations  (7a-7b):  the  latter  are  differentiated  with  respect  to  the  spatial  variable  in  this  effort. 
The  result  gives  for  the  space-charge  region  a  system  of  two  second  order,  ordinary,  non- 
linear, coupled  differential  equations  with  one  first  order,  coupled,  ordinary  differential 
equation.  The  non-dimensional  parameters  given  in  the  bracket  notation  are  listed  in  Tables 
2  and  10: 
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Essentially,  the  conservation  laws  are  combined:  species  momentum  now  reflects  the 
influence  of  charge  production  terms.  Here,  the  left-hand  side  of  equations  (10a- 10b)  is  the 
rate  of  change  of  the  number  density  gradient,  an  accelerating  behavior.  The  gradient  of  the 
species  number  density  is  found  in  the  first  term  on  the  right-hand  side  of  set  (10a- 10b),  a  rate 
influenced  by  the  strength  of  the  electric  field  as  measured  by  the  imbalance  between  the 
species  number  densities,  the  Gauss  equation  (10c).  The  first  term  of  the  right-hand  side  of 
(10a- 10b)  describes  the  interaction  of  coulombic  and  gradient  forces,  the  influence  of  the 
equation  of  motion  as  described  in  the  earlier  formulation  of  (1-5);  the  third  term  completes 
the  formulation  of  motion  in  this  model.  The  general  concept  of  energy  conservation  begins 
to  take  shape  with  the  second  and  fourth  terms  of  (10a- 10b):  The  former  has  included  two- 
body  recombination  effects  as  diminished  not  only  by  coulombic  interactions,  but  also  by 
diffusion.  The  fourth  term  illustrates  a  source  of  energy  to  the  system ,  ionization,  again  as 
influenced  by  diffusion.  Essentially,  the  accelerating  terms  on  the  left-hand  side  of  (10a)  and 
(10b)  are  described  by  balances  of  motion  (gradients)  and  energy  (ionization,  recombination) 
where  the  exchange  of  energies  is  brought  about  and  influenced  by  diffusion  processes.  It 
should  stand  to  reason  that  as  the  species  density  gradient  increases,  energy  losses  should 
become  more  dominant  until  the  gradients  diminish  their  contribution  to  the  flow: 
recombination  should  begin  to  take  hold  until  eventually  equilibrium  is  restored,  balancing 
diffusion  with  energy  source  and  sink. 
B.         TRANSITION  ZONE:  QUASI-NEUTRALITY 

In  order  to  determine  the  behavior  of  the  anode  fall,  boundary  conditions  have  to  be 
applied  to  equations  (7-9)  or  to  set  (10).  In  the  case  of  the  positive  column  and  the  cathode 
fall  region,  it  is  known  [Refs.  45,  46,  47,  48]  that  for  the  isothermal  case,  the  spatial 
derivatives  of  the  particle  density  and  the  drift  velocity  become  singular  if  the  drift  velocity 
attains  the  ambipolar  sound  speed  [Ref  49,  pp.  77-86;  50].  This  apparent  singularity  appears 
near  the  wall,  marking  the  end  of  the  quasi-neutral  region  and  the  beginning  of  the  space- 
charge  sheath  adjacent  to  the  wall;  the  ambipolar  solution  breaks  down  because  the  space 
charge  density  cannot  become  larger  than  the  local  density  of  ions  or  electrons.  [Ref.  50;  33; 
52:  pp    202-204].  Hence,  researchers  have  always  chosen  the  ion  sound  speed  as  the 
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appropriate  bound  for  the  ion  velocity.  This  condition  seems  a  natural  compromise  between 
Bohm's  criterion  for  the  collisionless  sheath  which  requires  the  ion  velocity  exceeds  or  is  at 
least  equal  to  the  ion  sound  speed  and  Persson's  criterion  for  an  ambipolar  flow  at  the  plasma 
boundary,  which  states  that  the  ion  velocity  is  in  fact  less  than  the  ion  sound  speed.  However, 
these  criteria  are  incompatible  as  was  shown  by  Godyak  and  Sternberg.  [Ref.  34]  For 
example,  in  the  numerical  treatment  of  the  plasma-wall  problem  without  plasma-sheath 
separation,  the  ion  velocity  reaches  ion  sound  speed  in  the  region  where  plasma  neutrality  is 
violated  [Ref.  51:  p.  51]. 

It  seems  that  since  there  is  no  strict  boundary  between  the  plasma  and  the  sheath,  it 
is  not  obvious  which  boundary  conditions  work  the  best;  much  controversy  remains  for  the 
negative  electrode.  Attempts  have  been  made  to  bridge  or  to  remove  these  apparent 
discontinuities  by  introducing  a  transition  layer  [Ref.  27]  or  by  choosing  various  plasma- 
sheath  boundary  conditions  [Refs.  45,  34].  For  example,  Godyak  and  Sternberg  [Ref.  34] 
show  that  the  usual  choice  of  choosing  the  ion  sound  speed  as  the  boundary,  leads  to 
discontinuities  of  the  potential,  plasma  density  and  velocity  gradients.  These  workers 
therefore  introduced  as  their  boundary  condition  a  specific  value  of  the  electric  field  that 
provides  a  breakdown  of  the  neutrality  at  the  plasma  boundary  and  so  treated  the  sheath  and 
ambipolar  model  separately,  attaining  a  smooth  transition  between  sheath  and  plasma  -  for  the 
cathode  region.  In  their  model,  however,  only  constant  ionization  was  assumed  for  species 
production. 

Fundamentally,  the  issue  of  transition  in  the  cathode  region  remains  somewhat 
controversial,  a  subject  even  more  unclear  when  the  positive  electrode  is  considered,  where 
essentially  scant  information  is  available.  But  a  transition  to  ambipolar  diffusion  must  occur 
[Ref.  52:  pp.  64-68],  a  transition  where  the  electric  field  domination  shifts  to  that  of  a  density 
gradient  driven  region,  an  area  where  now  the  length  scale  for  the  number  density  becomes 
important. 

To  capture  this  effect  and  to  retain  the  role  of  charge  production  in  the  dynamical 
behavior  of  the  system,  equations  (7)  are  again  differentiated  allowing  the  species  densities 
to  be  written  in  terms  of  the  production  terms.    Applying  the  quasi-neutral  condition, 
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n~ne    ,  then  gives  the  transition  model,  equation  (1 1)  which  can  also  be  obtained  from  the 


sheath  formulation,  set  (10). 
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Equation  (11) 

In  comparison  to  equations  (10a- 10b),  the  space-charged  model,  the  effects  of 
diffusion  are  clearly  seen  not  only  on  the  rate  of  change  of  the  number  density  or  density 
gradient  term  (first  term  of  (11)),  but  also  on  the  energy  balance  terms,  recombination 
(second  term  of  (1 1))  and  ionization  (third  term  of  1 1).  Furthermore,  coulombic  interactions 
are  no  longer  evident  in  the  energy  balance,  playing  only  a  secondary  role  in  the  density 
gradient  description  (first  term  of  (1 1)). 

A  key  character  of  equation  (1 1)  is  that  ambipolar  diffusion  is  not  yet  fully  formed, 
a  consequence  of  the  transition  region  in  which  perhaps  the  instability  of  the  sheath, 
electrostatic  repulsion  of  ions,  is  carried  forth  by  gradient  effects  so  that  ions  are  further 
pushed  from  the  positive  plate.  The  growth  in  densities  away  from  the  plate  then  affects  their 
free  diffusion,  a  process  evolving  into  ambipolar  diffusion. 
C.         AMBIPOLAR  REGION 

When  the  density  of  electrons  and  ions  becomes  large  enough,  their  mutual  Coulomb 
fields  affect  their  free  diffusion.  This  modification  can  be  written  in  terms  of  the  space  charge 
field  as  given  by: 


r       =  n,     v.     -±u.     n.     E-D.     Vn. 

ife  i(e     i,e  ^i,e     i,e  i,e         i,e 


Equation  (6) 


From  (6),  it  is  seen  that  should  the  species  flux  densities  not  be  equal,  a  significant  charge 
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imbalance  would  arise.  For  example,  if  the  plasma  is  much  larger  than  a  Debye  length,  it  must 
be  quasi-neutral  and  the  diffusion  rates  of  the  species  would  adjust  themselves  so  that  the 
same  rate  from  the  plasma  into  the  region  of  quasi-neutrality  is  achieved:  the  electrons  have 
a  larger  thermal  velocity  as  compared  to  the  ions  and  tend  to  leave  the  plasma  first.  A 
positive  charge  is  thus  left  behind  creating  an  electric  field  of  such  polarity  as  to  retard  the 
loss  of  electrons  and  to  accelerate  the  loss  of  ions:  the  plasma,  by  its  very  definition,  is  a 
neutral  fluid.  [Ref  38:  p.  139]  From  this  formulation,  the  ambipolar  diffusion  coefficient  can 

D  \L+Dp. 
be  determined  as       D  =  — — — —        where  the  Einstein  relation  may  be  applied  to  write 

mobilities  in  terms  of  diffusions.  It  should  be  pointed  out  that  the  ratio  D/u  represents  a 
measure  of  energy  [Ref.  52:  pp.  60-63],  and  that  if  Te=T; ,  as  is  assumed  in  this  work,  the 
ambipolar  electric  field  tends  to  enhance  the  diffusion  of  ions  by  a  factor  of  two,  a  diffusion 
that  is  primarily  controlled  by  the  slower  species:  Da~2Dj.  The  primary  function  of  the 
ambipolar  field  is  then  to  serve  as  an  energy-exchange  mechanism  for  transferring  random 
kinetic  energy  from  the  electrons  to  the  ions  such  that  the  ion  diffusion  continues  away  from 
the  positive  plate.  [Ref.  45] 

The  model  for  the  ambipolar  region  is  thus  formulated,  again  following  naturally  from 
equation  (10)  or  from  equations  (7-9),  the  sheath  model:  proceeding  as  before,  equations  (7) 
are  differentiated  with  respect  to  the  spatial  variable  allowing  set  (9)  to  contribute  to  the 
densities.  However,  unlike  previous,  now  set  (7)  is  multiplied  by  the  respective  species 
diffusion  coefficient,  the  result  which  gives  one  second  order,  ordinary  differential  equation 
for  the  density  in  terms  of  not  only  the  necessary  chemistry,  but  also  in  terms  of  the  ambipolar 
diffusion,  Da: 


2d2(n) 
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dp        LD;    De 
In  comparison  with  equation  (1 1),  it  is  immediately  seen  that  the  density  gradient 
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contribution,  the  rate  of  change  of  number  density,  is  no  longer  evident.  The  only  terms 
remaining  in  (12)  are  those  of  charge  production;  the  system  chemistry.  As  the  flow 
approaches  equilibrium  near  the  plasma  boundary,  the  number  density  must  approach  that  of 
Saha  so  that  the  definition  of  "plasma"  is  upheld. 

The  physical  model  in  summary:  the  electric  field  pushes  the  flow  from  the  electrode 
surface  toward  the  density  gradient  dominated  region  in  which  the  density  gradient  pushes  the 
system  toward  the  plasma  boundary  where  eventually  energy  conservation  must  be  satisfied. 
Ionization  must  balance  recombination  such  that  Saha  equilibrium  is  attained  at  the  plasma 
boundary.  It  is  in  this  region  that  the  largest  of  the  length  scales  has  influence:  the  current 
or  ionization  scale.  This  last  scale  defines  the  extent  of  the  anode  fall  region,  a  region  where 
flow  chemistry  becomes  paramount,  determining  plasma  equilibrium. 

But  is  the  anode  fall  stable?  Previous  attempts  [Ref  53]  at  numerical  analysis  have 
proven  challenging  at  best.  So  a  better  understanding  of  the  nature  of  equations  (7-9)  and 
sets  (10-12)  is  needed:  in  the  next  chapter,  non-linear  analysis  techniques  are  applied  to 
mathematical  models  presented  previously. 
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V.  ANALYSIS 

A.    LOCAL  ANALYSIS 

To  begin  the  analysis  of  the  anode  fall,  non-linear  methods  as  briefly  outlined  in 
Appendix  C  are  used:  further  information  on  particular  ideas  can  be  found  in  references  54 
and  55,  among  others. 

1.         Space-Charge  Region:  Sheath 

To  begin  an  analysis  of  either  equations  (7-9)  or  equations  (10),  non-linear  methods 
as  outlined  in  Appendix  C  are  used:  first,  equations  (10)  are  written  as  a  system  of  first  order 
differential  equations,  then  stagnation  points  within  both  systems,  equations  (7-9)  and  the  first 
order  formulation  of  equations  (10),  are  determined.  Next,  the  local  behavior  about  these 
stagnation  regions  in  terms  of  system  characteristics  is  evaluated,  leading  to  a  qualitative 
stability  picture  of  the  system  dynamics. 

a.  Equilibrium  (Stagnation)  Regions 

First  equations  (7-9)  and  equations  (10),  the  anode  sheath  models  are 
analyzed  (non-dimensional  parameters  are  given  in  Table  2): 
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The  original  model,  equations  (7-9)  can  immediately  analyzed,  whereas  equations  (10) 
have  to  be  recast  as  a  set  of  coupled,  first-order  differential  equations. 
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Equation  (10a) 

Equation  (10b) 
Equation  (10c) 


To  find  the  stagnation  points  within  these  first-order  system,  the  derivative  terms  are  each  set 
to  zero  and  non-dimensional  parameters  from  Table  2  are  applied.  The  results  give  for  the 
equilibrium  criteria: 

1 .   From  equations  (8)  and  (10c),  ni*=ne*=n*     as  expected;  this  result  is  then 
incorporated  into  equations  (9a-9b)  as  well  as  (10a- 10b). 


v  v 


2.  Similarly,  from  equations  (9a-9b)  and  (10a- 10b),      n  *  -  0      or    n  *  =  — - —    ;  a 


2     o 


stagnation  point  appears  when    n*  =n*  -n*  =  0      and  when  the  plasma  boundary 

is  approached.  The  latter  condition  immediately  emphasizes  the  flow  chemistry, 
ionization  and  recombination.  Without  recombination,  the  second  stagnation  would 
not  be  present,  as  seen  from  either  equations  (9a-9b)  or  equations  (10a- 10b)  with 
d/dy[  ]=0.  The  former  equilibrium,  on  the  other  hand,  seems  artificial  at  first,  unless 

the  species  densities  behave  such  that  they  cross  when      n*=n*-nm=0 
immediately  implying  a  quasi-neutral  condition  since  equation  (8)  or  equation  (10c) 

require    n.*  -ne* =n*    at  stagnation.  In  general,  the  fixed  points  found  stem  purely 

from  the  chemistry  behavior,  equations  (9a-9b)  and  consequently  (10a- 10b)  in  which 
ionization  and  two-body  recombination  are  the  only  processes  included.  If,  on  the 
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other  hand,  three-body  recombination  is  added  to  the  equations,  then  additional 
stagnation  regions  become  evident  since  equations  (9a-9b)  are  now  cubic.    The 

equilibrium  region  near      n*  -n*  =n*  =  0     still  remains,  however.    Similarly,  if 

ionization  is  the  only  production  term  considered  in  the  sheath  model,  then  the 
stagnation   region   near  the   plasma   boundary  "disappears."      The   condition 

n*  =n*  -n  * = 0     remains  for  all  cases:  ionization  only,  ionization  with  two-body 


recombination,  ionization  with  two-  and  three-body  recombination.  So  it  appears  that 
continuity  (equations  (9))  requires  flow  equilibrium  when  perhaps  quasi-neutral 
conditions  are  initially  met,  at  the  sheath  boundary,  and  when  the  plasma  interface  is 
approached.  It  is  recombination  that  establishes  the  latter  while  the  presence  of 
ionization  stipulates  the  former. 


3.  The  current  densities  at  the  equilibrium  position  can  be  found  by  using    n  *  =  0 
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with    equations    (7a- 7b):     if         n*  =  — - —      ,    then 
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since     n   -n    -n         when    n,  =«,  =n 


Here,  the  influence  of  not  only 

ionization  and  recombination,  but  also  of  diffusion  on  the  current  density  distributions 
are  clearly  seen.  Similarly,  from  equations  (7a-7b),    jie    =0  V«*=0 


4.  Charge  conservation  is  used  to  isolate  the  value  of  the  electric  field  at  stagnation 

vov 
near      the      plasma      boundary,      where        w*  = .  That      is, 
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from  which   E  *  is  obtained  graphically  since  the  right-hand  side  of  this  expression 
simply  represents  a  constant.  Applying  the  Arrhenius  function  [Ref.  52:  pp.  94-109] 


for  the   ionization,     v~e 


=  —        and   using   numerical   data  for   nitrogen 


yields  £**0.54   .  From  this,     j/  *  0.99816,  j*»  0.0005445,  »*«0. 15695       . 

These  results  show  the  relative  magnitude  of  the  electron  and  ion  currents  within  the 
anode  fall,  je»jj  as  expected.  If,  on  the  other  hand,    j, ,e *  =0  V n  *  =  0    is  used,  an 

inconclusive  result  is  obtained:     E  *    cannot  be  isolated  from  either  equations  (7-9) 


or  equations  (10).  Further,     v  ~  e L    ->=  —     does  not  permit  zero  field.  Therefore, 


it  seems  reasonable  that  a  finite  field  should  exist  over  the  entire  range  of  the  anode 
sheath.  Von  Engel  [Ref.  56:  pp.  12-20]  highlights  such  arguments:  a  finite  field  must 
be  allowed,  which  could  then  introduce  a  potential  trough.  A  virtual  anode  is  thought 
to  exist  if  the  space  charges  have  the  same  sign  as  the  nearby  electrode,  with  the 
emitting  surface  located  at  the  potential  minimum  [Ref.  56:  pp.  14-17]. 

The  fixed  points  found  in  the  model  represented  by  equations  (7-9)  and  in  part  by  set  (10)  are 
listed  in  Table  5  below,  featuring  the  role  of  ionization,  recombination  and  diffusion  in  the 
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existence  of  the  fixed  point  near  the  plasma  boundary.  Again,  continuity  (equations  (9)) 
seems  to  require  flow  equilibrium  when  perhaps  quasi-neutral  conditions  are  initially  met,  the 
sheath  boundary,  and  when  the  plasma  interface  is  reached.  It  is  the  recombination 
phenomena  that  establishes  the  latter  stagnation  region  while  the  presence  of  ionization 
stipulates  the  former. 


Variable 

V 

Near  the 

Sheath 

Boundary 

»2* 

Near  the  Plasma  Interface 

Numerical 
Value  for 
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«' 
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a2"o 
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Table  5:  Stagnation  Regions  in  General  Sheath  Model  -  Ionization  with  Two-body 

Recombination 
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b.         Jacobian  at  the  Equilibrium  Positions 

To  determine  the  system  characteristics  about  its  fixed  points,  a  local  Taylor 
series  expansion  is  carried  out  about  the  equilibrium  of  interest  (Appendix  C);  the  rate  of 
change  of  the  flow  is  given  in  terms  of  the  Jacobian,  which  for  equations  (7-9)  becomes: 
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Here,  six  variables  affect  the  system  dynamics:  ion  and  electron  number  densities, 
electric  field,  species  difiusivities,  and  ionization.  Of  these,  the  difiusion  coefficients  are 
assumed  constants,  so  essentially  four  variables  drive  the  behavior:  species  densities,  electric 

field  and  consequently  ionization.   If     Ainner     is  evaluated  at  the  stagnation  points,  the 

resulting  eigenvalues  describe  the  dynamics  in  the  neighborhood  of  the  fixed  point  (Appendix 
C).  Specifically,  at  the  plasma  boundary  equilibrium  with  nitrogen  data: 


A,  =  1.376704 
X3  =  0.630455 
A5= -0.63279 


A2  =  0 

X4=  -1.37437 


For  this  case,  the  eigenvalues  are  real  and  distinct,  leading  to  a  multi-dimensional  type  of 
"saddle"  behavior  (Appendix  C).  The  stagnation  region  near  the  plasma  boundary  is  unstable: 
the  largest  eigenvalue  is  positive,  indicative  of  exponential  growth.  Appendix  D  lists  the  roots 
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of  the  characteristic  polynomial  of    Ainner     .  It  appears  Xx  gives  exponential  growth  -  the 

process  of  ion  repulsion  from  the  positive  plate  forces  the  instability  of  the  sheath.  The  sheath 
drives  the  system  toward  the  first  stagnation  criterion,  toward  quasi-neutrality. 

Provided  a  finite  field  exists  when     n*  -n*  =n  *  =  0     ,  then  the  fixed  point  in  its 


vicinity  gives  a  similar  qualitative  behavior:  the  eigenvalues  of  the  resulting  Jacobian  are  given 
by 

Xx=[a]E         ^3=0 
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These  eigenvalues  are  always  real  with  two  characteristics  equal  while  the  remaining  are 
distinct:  Table  6  lists  some  numerical  examples  for  nitrogen  data. 
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2                         2 
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-0.00053 
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([c]£J2-4 
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Table  6:  Eigenvalues  of  Stagnation  Region  ^=^=^1=0  with  Finite  Field  -  Numerical  Values 

for  Nitrogen 

Again,  the  essential  feature  of  the  sheath  is  instability:    ions  are  driven  from  the 
positive  wall  as  required  by  electrostatics.  The  stagnation  region  near     n*  -n*  -n*  =  0     is 

unstable,  perhaps  pushing  the  system  away  toward  the  recombination  induced  equilibrium 
near  the  plasma  boundary.  It  is  the  strength  of  the  initially  applied  electric  field  that  imposes 

the  stability  characteristics  for  the  fixed  point     n* -n*  =n*  =0         This  leads  to  the 

following  question:  is  the  sheath  model,  as  presented  by  either  equations  (7-9)  or  set  (10), 
an  example  of  a  system  in  which  a  certain  type  of  "fixed  point"  is  always  present,  a  stagnation 
region  satisfying  quasi-neutrality  and  thereby  representing  the  sheath  boundary?  If  so,  does 
this  equilibrium  present  a  type  of  transcritical  bifurcation  (see  Appendix  C)  whose  stability 
perhaps  depends  on  the  initially  applied  electric  field?  The  foregoing  does  suggest  so. 

The  precarious  nature  of  the  fixed  points  found  in  the  sheath  model  suggest  that 
numerical  integration  methods  will  not  prove  fruitful  over  the  entire  domain  of  the  anode  fall. 
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However,  by  exploring  the  various  regimes  of  equations  (7-9)  in  the  form  of  the  second  order 
differential  equations  given  by  equations  (10),  further  insight  into  the  anode  fall  could  be 
expected. 

2.  Transition  Zone:  Quasi-Neutrality 

Previously,  it  has  been  shown  that  the  sheath  is  inherently  unstable,  that  the  ions  are 
electrostatically  repelled  from  the  electrode  (anode)  surface.  If  through  that  process  quasi- 
neutrality  occurs,  then  for  steady-state  conditions  either  density  gradients  or  diffusive 
processes  have  to  drive  the  system  toward  the  plasma  boundary  since  the  condition  of 

n.~n      stipulates  the  electric  field  be  constant;  so  something  else  has  to  continue  driving 

the  ions  away  from  the  anode  wall.  Equation  (11)  highlights  the  density  gradient  and  is 
considered  next  as  a  set  of  coupled  first-order  differential  equations: 
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Equation  (11) 


where  the  electric  field  is  constant  due  to  quasi-neutral  conditions,     ni ~ne 
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a.  Equilibrium  (Stagnation)  Regions 

To  find  the  stagnation  points  within  these  first-order  system,  the  derivative 
terms  are  each  set  to  zero  and  non-dimensional  parameters  from  Table  2  are  applied.  The 
results  gives  for  the  equilibrium  criteria: 

1.  Setting  d  [  ]/dy=0  in  equation  (11)  immediately  gives  the  previously  found 

v  v 

stagnation  condition  near  the  plasma  boundary:       n*  = The  fixed 

2     o 

point     n  *  =  0     again  seems  to  present  a  mathematical  phenomenon  since  equation 

(11)  does  not  describe  the  sheath  condition.   If  quasi-neutral  arguments  are  again 

used,  then  perhaps  the  condition    n  *  =  0     may  have  meaning  -  the  plasma  shields 

itself  against  external  influences  (electric  field),  and  so,  quasi-neutral  characteristics 

are  fundamental  to  its  definition.  The  stagnation  region  at     n  *  =  0    does,  perhaps, 

present  a  physical  phenomenon:  the  onset  of  quasi-neutrality  or  the  sheath  boundary, 
found  as  a  consequence  of  ionization  present  in  the  sheath  model  and  in  equations 
(11),  the  "transition"  model. 

2.  The  magnitude  of  the  constant  electric  field,  and  therefore  constant  ionization,  at 
the  start  of  the  quasi-neutral  transition  is  determined  by  the  matching  conditions  when 

n~ne    .  These  criteria  are  either  found  through  analytically  matching  equations  ( 1 0) 

with  equations  (1 1)  for  some  small  parameter(s),  or  numerically.  Numeric  methods 
will  be  considered  later  in  this  work.  Suggestions  toward  an  analytic  matching 
attempt  are  provided  in  paragraph  B  below. 

b.  Jacobian  at  the  Equilibrium  Positions 

To  determine  the  system  characteristics  about  its  fixed  points,  a  local  Taylor 
series  expansion  is  carried  out  about  the  equilibrium  of  interest  (Appendix  C);  the  rate  of 
change  of  the  flow  is  given  in  terms  of  the  Jacobian,  which  for  equation  (11)  becomes: 
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Since  the  electric  field  becomes  constant  at  the  onset  of  quasi-neutrality  at  the  sheath 
boundary,  the  electric  field  thus  locks  in  the  behavior  of  the  transition  region.  To  determine 

these  stability  characteristics,  the  general  eigenvalues  of  A,       .,.       are  determined: 
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where  the  non-dimensional  parameters  are  listed  in  Tables  2  and  10.    These  eigenvalues 
determine  the  dynamical  behavior  of  the  transition  when      A,       .,.         is  evaluated  at  a 

3  transition 

fixed  points,  albeit  the  electric  field  strength  is  determined  from  the  sheath  as  the  system 
reaches    n~ne    .  By  inspection  of  the  first  term  found  in  these  eigenvalues,  the  following 


is  observed: 

1 .  If  the  field  is  positive,  then  the  system  displays  exponential  decay,  stability;  but 
electrostatics  demands  that  ions  be  repelled  from  the  positive  plate.  This  leads  to  the 
expectation  that  the  transition  zone  is  also  unstable,  similar  to  the  sheath. 
Consequently,  the  field  should  be  negative  at  the  point  of  quasi-neutrality. 


2.     Using  the  stagnation  criteria  for      n  *  =  0     with  finite  field  indicates  that  the 
characteristics  should  display  stable  behavior  for  positive  field,  the  diffusion  difference 
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(DrDe)  is  negative,  and  unstable  behavior  for  negative  field.  The  initially  applied 
electric  field  strength  then  does  determine  the  stability  of  the  fixed  point  n  *  =  0  , 
an  equilibrium  which  appears  when  ionization  is  present  in  the  model. 

v  v 

3 .  To  find  the  behavior  near  the  plasma  boundary  when     n  *  = ,  nitrogen 

a2no 
data  along  with  an  arbitrary  selected  field  strength,    E(y-  -  )  =  -5.00   are  used  to 

i       e 

give  the  eigenvalues  Aj—0.041  and  A2=5.026:  the  eigenvalues  are  real,  distinct  and 
predict  a  saddle-type  behavior  near  the  plasma  boundary  -  a  robust  linearization  (see 
Appendix  C).  For  this  case,  the  flow  field,  locally,  is  characterized  by  a  very  strong 
unstable  manifold  corresponding  to  A2.  A  similar  result  was  found  by  analyzing  the 
sheath  model  presented  by  either  equations  (7-9)  or  equations  (10). 

4.  On  the  other  hand,  if  the  sheath  behavior  at  quasi-neutral  formation  gives  a 
positive  field,  e.g.,    E{y-  _-  )  =  5.00  ,  then  the  Jacobian  eigenvalues  predict  stability 

locally  about  the  fixed  point:  ^=0.021  and  Xj=-5.013  where  now  the  saddle  behavior 
is  marked  with  a  strong,  stable  manifold,  contrary  to  the  sheath  analysis  earlier. 

Thus,  the  criteria  found  at  the  onset  of  quasi-neutrality  at  the  sheath  boundary  form 
the  initial  conditions  imposed  upon  the  transition  zone  and  hence  govern  the  stability 

_  _       v  v 

characteristics  of  the  fixed  points  at      n*  =  0    and     n*  =  — - — 

aino 


3.         Ambipolar  Region 

In  the  sheath,  ions  are  repelled  electrostatically  until  quasi-neutrality  occurs,  where 
chemical  processes  balance  such  that  the  ion  density  approaches  that  of  the  electrons  -  the 
electric  field  becomes  constant.  If  this  balance  does  not  reflect  total  charge  neutrality,  i.e., 
Saha  equilibrium,  then  other  processes  (the  field  is  now  constant)  must  drive  the  system 
toward  such  equilibrium;  however,  a  stabilizing  mechanism  must  be  available  once  the 
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densities  are  near  the  plasma  boundary.  So  as  the  species  density  gradient  increases,  energy 
losses  should  become  more  dominant  until  the  gradients  diminish  their  contribution  to  the 
flow:  recombination  should  begin  to  take  hold  until  eventually  equilibrium  is  restored, 
balancing  diffusion  with  energy  source  and  sink  -  ambipolar  diffusion  becomes  the  vehicle 
which  takes  the  flow  density  toward  the  plasma  stagnation.    To  translate  these  physical 

processes  to  a  mathematical  model,  quasi-neutral  conditions,    n~ne    ,  are  applied  to  the 

sheath  equations  (10)  from  which  equation  (12)  results  naturally,  highlighting  the  ambipolar 
diffusion  in  its  full  form: 
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where  the  electric  field  is  constant  due  to  quasi-neutral  conditions,   n~ne    . 

a.  Equilibrium  (Stagnation)  Regions 

To  find  the  stagnation  points  within  this  first-order  system,  the  derivative 
terms  are  each  set  to  zero  and  non-dimensional  parameters  from  Table  2  are  applied.  The 
results  gives  for  the  equilibrium  criteria: 


1.    Setting  d  [  ]/dy=0  in  (12)  immediately  gives  the  previously  found  stagnation 


condition  near  the  plasma  boundary:      n  * 


VoV 


2     o 
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2.      Similar  to  the  transition  zone,  the  fixed  point   n*  =  0   again  presents  a 
mathematical    phenomenon   until   physical    arguments   are   applied   as   before. 
n  *  -  0    perhaps  indicates  the  onset  of  quasi-neutral  conditions. 

The  magnitude  of  the  constant  electric  field,  and  therefore  constant  ionization,  at  the 
start  of  the  quasi-neutral  ambipolar  region  continues  from  the  transition  zone  and  again  locks 
in  the  behavior  of  the  ambipolar  region:  the  initial  field  strength  determines  the  initial 
conditions  for,  and  consequently  the  stability  of,  the  transition  zone.  The  way  in  which  the 
density  gradient  then  develops  influences  the  initial  conditions  to  the  ambipolar  region.  All 
in  all,  the  initially  applied  electric  field  strength  impacts  the  stability  of  the  anode  fall,  in 
entirety  -  the  electric  field  directly  influences  ionization  and  hence  the  stagnation  found  near 
the  onset  of  quasi-neutrality. 

b.  Jacobian  at  the  Equilibrium  -  Analysis 

To  determine  the  system  characteristics  about  its  fixed  points,  a  local  Taylor 
series  expansion  is  carried  out  about  the  equilibrium  of  interest  (Appendix  C);  the  rate  of 
change  of  the  flow  is  given  in  terms  of  the  Jacobian,  which  for  equation  (12)  becomes: 
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If  the  stagnation  conditions  near  the  plasma  boundary  are  applied,  it  can  be  easily  be  seen  that 
the  ambipolar  eigenvalues  are  always  real  and  distinct,  again  predicting  the  saddle-type 
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behavior  (Appendix  C).  In  comparison  to  both  the  sheath  and  the  transition  zone,  the  electric 
field  no  longer  is  found  explicitly  part  of  the  Jacobian  eigenvalues,  however  an  implicit  role 
is  still  evident  in  terms  of  the  ionization.  Furthermore,  diffusion  and  chemistry  are  part  in  all 
eigenvalues,  from  sheath  to  ambipolar. 

Proceeding  as  before,  the  Jacobian  is  evaluated  at  the  plasma  boundary  stagnation, 

_      v  v 

n  *  = — - —    ,  using  nitrogen  data.  The  initial  field  strength  is  chosen  arbitrarily;  "initial"  is 

aino 

defined  as  that  condition  found  at  the  onset  of  the  ambipolar  region  as  determined  by 
matching  criteria  between  sheath  and  transition,  at  the  beginning  of  quasi-neutrality.  The 
resulting  Jacobian  eigenvalues  or  characteristics  are  again  robust  (system  non-linearities  do 
not  change  the   qualitative  nature  of  the  flow  )  [Ref.   54:   p.    163].      That  is,  if 

E(y-_-  )  =  -5.00   ,  Au=±  37.13,  whereas  if  E(y-_-  )  =5.00  ,  Xu  =  ±  30.4.  Although 

the  eigenvalues  are  real  and  distinct,  they  are  of  the  same  magnitude:  exponential  growth  is 
exactly  offset  by  exponential  decay  -  a  hallmark  of  energy  conservation  or  homoclinic 

behavior  (see  Appendix  C).  The  phase  portrait  for  equations  (12),      n(y)  vs  _         , 

dy 

using  nitrogen  data,  is  shown  in  Figure  (4). 
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Figure  4:  Phase  Portrait  for  Ambipolar  Region,  Equation  (12),  with  Nitrogen  Data  for  64 
Different  Initial  Conditions  using  Ionization  and  Two-Body  Recombination 

Can  the  stagnation  near  the  plasma  boundary  be  reached?  Is  this  region  stable?  The 
behavior  of  the  ambipolar  region  is  driven  by  the  initial  conditions  as  applied  to  equations  (12) 
as  shown  by  figure  (4):  if  the  initial  conditions  are  such  that  region  I  is  dominant,  then  the 
manifolds  of  the  stagnation  near  the  plasma  boundary,  the  saddle  behavior  (see  Appendix  C), 
confine  the  system  trajectories.  If,  on  the  other  hand,  the  initial  conditions  are  such  that 
region  II  is  dominant,  then  the  system  will  become  unstable.  The  manifolds  of  this 
linearization  form  the  threshold,  at  least  locally,  which  perhaps  determines  the  requirements 
for  quasi-neutrality  -  the  homoclinic  trajectories,  as  given  below: 
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Here,  the  contributions  of  ambipolar  diffusion,  ionization  and  two-body  recombination  are 
clearly  seen;  two-body  recombination  is  found  in  the  non-dimensional  parameters  [h]  and  [j]. 
More  often  than  not,  these  homoclinic  trajectories  can  be  volatile:  any  numeric  round-off 
error  incurred  through  whatever  computational  algorithm  used,  will  induce  sufficient  non- 
linearities  that  will  prevent  the  solution  from  being  "pinned"  to  these  manifolds,  the  manifolds 
will  break.  Even  if  it  were  possible  to  achieve  this  exact  trajectory,  the  initial  conditions 
imposed  by  matching  the  transition  region  to  the  ambipolar  equation  will  most  likely  make 
such  efforts  challenging.  That  is,  the  initial  conditions  provided  to  the  quasi-neutral  region 
by  the  sheath  determine  the  behavior  of  the  ambipolar  evolution:  either  stable  periodicity  as 
shown  by  region  I  in  Figure  (4),  or  exponential  growth  as  indicated  by  region  n.  In  either 
case,  the  behavior  of  the  sheath  ultimately  determines  the  behavior  of  the  outer  regions  and 
therefore  the  anode  fall,  in  general. 

4.         Comments/Observations 

The  preceding  suggests  that: 

1 .  Ions,  created  through  electron  bombardment  of  neutral  atoms  throughout  the 
anode  fall  region,  are  repelled  electrostatically  from  the  positive  electrode  toward  a 
location  of  quasi-neutrality  at  the  sheath  boundary  as  well  as  toward  an  unstable 
stagnation  region  near  the  plasma  interface.  These  equilibria  form  as  a  result  of  the 
chemistry  present.  That  is,  continuity  (equations  (9))  seems  to  require  flow 
equilibrium  when  quasi-neutral  conditions  are  initially  met  within  the  sheath,  and  when 
the  plasma  boundary  is  reached.  It  is  recombination  that  establishes  the  latter  while 
the  presence  of  ionization  stipulates  the  former.  Fundamentally,  the  sheath  as 
formulated  is  inherently  unstable  -  ions  are  driven  out  of  the  system,  toward  the 
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negative  electrode:  it  is  the  role  of  the  anode,  after  all,  to  provide  the  positive  current 
to  the  cathode  [Ref.  43:  p.  59]. 

2.  At  the  moment  of  quasi-neutrality,  when    n~ne    ,  there  appears  a  conservative 

phenomenon,  such  that  the  energy  growth  is  matched  by  the  energy  loss  -  a  chemical 
balance.  Is  the  sheath  model  as  presented,  an  example  of  a  system  in  which  a  certain 
type  of  "fixed  point"  is  always  present,  a  stagnation  region  satisfying  quasi-neutrality? 
If  so,  does  this  equilibrium  present  a  type  of  transcritical  bifurcation  whose  stability 
depends  on  the  initially  applied  electric  field?  The  brief  analysis  given  in  this  chapter 
does  suggest  such  behavior  -  the  onset  of  quasi-neutral  characteristics,  the  sheath 
boundary,  might  be  indicative  of  a  transcritical  bifurcation. 

3.  The  quasi-neutral  balance,  however,  is  insufficient  to  provide  the  necessary 
stability  condition  to  ensure  the  system  progresses  toward,  and  remains  in  the  close 
vicinity  o£  the  plasma  boundary  fixed  point.  As  such,  any  numeric  solutions  will  most 
likely  be  difficult  to  achieve.  That  is,  the  initial  conditions  provided  to  the  quasi- 
neutral  region  by  the  sheath  determine  the  behavior  of  the  ambipolar  evolution:  either 
stable  periodicity  as  shown  by  region  I  in  Figure  (4),  or  exponential  growth  as 
indicated  by  region  n.  In  either  case,  the  behavior  of  the  sheath  ultimately  determines 
the  behavior  of  the  outer  regions  and  therefore  the  anode  fall,  in  general. 

B.         ANALYTICAL  MATCHING  -  SUGGESTIONS 

The  analytic  analysis  presented  in  this  chapter  suggests  that  the  anode  fall  model  used 
in  this  work  does  not  represent  sufficient  dissipation,  preventing  the  system  from  adjusting 
itsel£  to  adhere,  to  the  "stable  manifold"  of  the  stagnation  region  near  the  plasma  boundary. 
The  model  is  most  likely  incomplete  in  its  physical  representation.  Perhaps  analytically 
matching  the  second  order  sheath  equations  to  the  transition  set  which  are  then  matched  to 
the  ambipolar  equation,  may  reveal  an  inconsistency  between  some  of  the  physically 
represented  terms.  Albeit  this  process  seems  challenging  at  first,  perhaps  some  of  the  ideas 
suggested  next  could  be  applied: 

1.  Three  length  scales  have  been  found  when  equations  (7-9)  were  non- 
dimensionalized  (Appendix  A  and  Table  3).  Three  "small"  parameters  could  then  be 
used  to  match  equations  (10-12)  as  follows:  Assume  the  sheath  equations  (10) 
behave  rapidly  over  a  very  small  distance,  hence  a  fast  changing  parameter  should  be 
applied,  1/5=0(1^).  Further,  assume  the  ambipolar  model  (12)  as  varying  slowly  in 
comparison,  e~0(Ln  or  Lg).  In  between  is  a  region  of  transition  where  r|=e/6.  But 


42 


where  are  these  equations  matched? 

3.  To  determine  the  location  of  matching,  the  fixed  points  found  from  equations  (7-9) 
could  be  used  as  a  starting  point.  The  linearization  of  equations  (10-12)  about  these 
stagnation  regions,  then  provides  the  criteria  for  matching:  the  respective  manifolds 
describe  the  number  density  and  the  density  gradients  evident  near  these  fixed  points. 
The  interpretation:  the  manifolds  given  by  the  linearization  about  the  fixed  points 
gives  the  conditions  for  matching,  criteria  which  must  be  met  by  the  preceding  region 
as  the  system  is  driven  into  the  next  zone.  The  manifolds  for  the  fixed  point  near  the 
plasma  boundary  have  already  been  presented  in  section  A3,  above.  The  manifolds 

for  the  fixed  point  near  the  sheath  boundary  can  be  found  from      ^transition      by 

using  any  symbolic  mathematic  utility  such  as  MAPLE  or  MATHCAD.  The 
following  is  an  illustration  of  the  above  concept: 


Transition  Model:  Am  bipolar  Model: 

Linearization  Manifolds  Linearization  Manifolds 

Determine  Sheath  to  m  ^  Determine  Transition  to 

Transition  Match     .  J^f^^  Ambipolar Match 

t7c  fS            StaS**ti°*  Boundary 

V                ^    am   **       4iT  Stagni 


g 


Stagnation  e„«,.,Wfl/ 


Distance  from  Wall 

3.  The  initial  conditions  that  are  used  for  the  sheath,  however,  remain  a  "shooting" 
problem.  Perhaps  a  numeric  iteration  scheme  where  many  initial  conditions  are  tested 
might  reveal  some  pattern  for  better  selecting  these  conditions. 
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VI.  NUMERICAL  RESULTS  FOR  NITROGEN 

Using  available  MATLAB  4-5-stage  Runge-Kutta  algorithms,  both  equations  (7-9) 
and  equations  (10),  as  first  order  differential  equations,  are  integrated  numerically.  The  data 
thus  obtained  were  then  exported  into  a  Quattro-Pro  spreadsheet  for  plotting  and  easier 
analysis. 

1 .  Length  scales  (Table  3)  to  determine  the  non-dimensional  parameters: 
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2.    Non-dimensionalized  initial  conditions  (Table  4  lists  the  numerical  data  for 
nitrogen): 
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The  drift  currents  are  used  to  generate  the  appropriate  initial  conditions  at  the 
electrode  wall  since  at  the  wall,  the  field  induced  motion  will  exceed  that  induced 
thermally;  it  is  assumed  that  initially  any  density  gradients  are  zero.  To  gain  insight 
as  to  the  impact  of  initial  conditions  on  the  system,  the  initial  electric  field  is  varied 
from  0.01  to  20(1 04)  V/m.  The  results  are  given  below.  Of  note,  the  method  used 
to  non-dimensionalize  the  variables  (see  Appendix  A)  gives  rather  small  "non- 
dimensional  numbers"  in  terms  of  magnitude,  especially  when  the  integration  results 
shown  are  "close"  to  the  electrode  wall. 
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A.         SPACE-CHARGE  REGION,  SHEATH 

In  overview:  For     E(0)  < Ethreshold  =  0.04     ,  the  ion  density  behaves  as  expected, 

increasing  monotonically,  being  primarily  influenced  by  its  gradient.  Both  the  low  initial  field 
strength  and  the  ion  density  gradient  are  insufficient,  however,  to  push  the  species  toward 

quasi-neutrality  when  nj-ng.   When     E(0)>  Ethreshold  =  0.04    ,  the  initial  field  becomes 

strong  enough  to  affect  the  ion  density  as  compared  to  its  gradient;  the  ions  are  pushed 
toward  what  appears  as  a  density  trough.     Furthermore,  for  certain  initial  values  of 

A»  =  \nf  -  ne\       vs    E(0)      ,  the  ions  are  pushed  past  their  minimum  density  point, 

eventually  approaching  the  electron  density  distribution,  after  which  both  densities  remain 
close  in  magnitude  for  a  short  distance;  essentially,  for  adequately  strong  initial  electric  field, 
quasi-neutrality  is  formed  for  a  brief  distance.  The  quasi-neutral  characteristic  is  inherent 
within  the  energy  balance  found  in  the  ions,  the  balance  between  ionization  and 
recombination;  without  this  balance,  the  space-charge  region  remains. 

1.  Varying  Initial  Conditions:  Impact  on  System  Behavior 

To  gauge  the  impact  of  various  initial  conditions,  specifically  combinations  of  species 
density  and  electric  field,  on  the  sheath  dynamics,  equations  (10)  were  integrated  for 
^(O^IO8  [1/m3],  ne(0)=1012  [1/m3]  and  E(0)=6(104)  [V/m]  using  the  length  scales  as  given 
in  Table  3.  These  conditions  were  chosen  after  varying  the  initial  conditions  as  will  be 
described  below. 

Specifically  for  the  ion  density  in  the  sheath: 
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Ion  Density  vs  Distance  from  Wall 

Varying  Initial  Electric  Field 


1.10E-11 


Ei 

T3 
i 

§1 


'55 

eg 

© 

Q 

c8 

o 


05E-1 1 
00E-11 
50E-12 
00E-12 
50E-12 


8.00E-12 


,=i 

.1. 

™T 

j 

■" 

For  sufficiently  strong 
initial  field,  ions  are  pushed 
toward  density  trough 

i 

0.00E+00    2.00E-07    4.00E-07    6.00E-07    8.00E-07    1.00E-06    1.20E-06 

Distance  from  Electrode  Wall  [m] 


E^O)  [non-dim]        „     3.Q2E  02    ■     4.00E-02    ne(°):  1QA12  [l/m3],  ni(0):  10A8  [l/m3] 


Figure  5:  Ion  Density  vs  Distance  -  E(0)  <  E^^y 

Figure  (5)  shows  the  ion  density  vs  distance  from  the  electrode  wall  for  distances  very 
close  to  the  electrode  surface  -  in  comparison,  see  Figure  (8)  given  later  in  this  chapter. 
Again,  the  non-dimensionalization  used  makes  the  magnitudes  for  the  density  rather  small, 
but  these  magnitudes  represent  non-dimensional  quantities  and  so  should  be  multiplied  by  n0 
(see  Table  4). 

At  the  electrode  wall,  in  its  immediate  vicinity,  there  are  at  least  two  major  vehicles 
which  drive  the  system  away  from  the  wall  in  steady-state:  the  initial  electric  field  strength  and 
the  density  gradient,  such  that  if  the  field  strength  is  low  enough,  the  gradient  dominates  the 
system  whereas  if  the  field  strength  is  high  enough,  the  gradient  influence  becomes  secondary 
in  nature.  Figures  (6)-(7)  depict  the  magnitude  of  the  ion  density  gradient  vs  distance  from 
the  electrode  wall  and  the  electric  field  vs  distance,  respectively.  Comparison  with  Figure  (5) 


leads  to  the  speculation  that  for  low  electric  field  strength,     E(0)<Ej 
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ion  density  is  governed  by  its  gradient,  whereas  when   E{  0 )  >  Ethre  .  u  =  0.04    ,  the  electric 


field  drives  the  ions  from  the  wall.  When    E(0)  =  0.04      the  boundary  is  reached  where 


initial  density  gradient  effects  are  surpassed  by  those  of  the  field  (Figures  (5)  and  (6)) .  As 
shown  in  Figure  (5),  the  ion  density  appears  to  increase  from  the  electrode  wall  for  low  initial 
field;  yet  this  growth  does  not  continue  for  very  long.  Essentially,  the  field  strength  is 
insufficient  to  drive  the  density  further  while  the  density  gradient  is  too  weak  to  provide  the 
necessary  impetus  for  the  ion  and  electron  densities  to  match.  In  general,  Figures  (5)-(7) 
show  results  for  integrated  variables  in  close  vicinity  of  the  electrode  surface  -  in  comparison, 
see  Figure  (8). 

Ion  Density  Gradient  vs  Distance 

Varying  Initial  Electric  Field 


.§  4.00E-14 
"O 

C  2.00E-14 
O 

So.ooe+oo 

C 

•2-2.0OE-14 
T3 

2-4.00E-14 
.&-6.00E-14 

m 

HJ-8.0OE-14 

Q 

C-1.00E-13 

O 


±-           .     ... 

i 

E(0)<=0.04  [non-dim] 

J^*^ 

•""i 

..  .   -       _.  _-_  -    . 

j 

I 

j         i 

p^v 

j 

For  sufficiently  strong  initial  field,  ions 
are  pushed  toward  density  trough.  If 
low  field,  gradient  drives  density. 

■ 

^vj 

l\ 

- 

!     I     i     I     I     I     I 

0.00E+00    2.00E-07    4.00E-07    6.00E-07    8.00E-07    1.00E-06    1.20E-06 

Distance  from  Electrode  Wall  [m] 


E$)  [non-dim]     __  3.92E-02—  4.00E-02    "e(°):  1QA12  [1/m3l.  ^)-  10Ag 


Figure  6:  Ion  Density  Gradient  vs  Distance  -  E(0)  <  E,^,^ 
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Electric  Field  vs  Distance 

Varying  Initial  Electric  Field 
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Figure  7:  Electric  Field  vs  Distance:  E(0)  <  E^g,^ 
A  density  minimum  occurs  whith  the  initial  electric  field  strength  increases  (Figure 
(8)).  Moreover,  the  density  decline  reaches  a  minimum  according  to  the  electric  field  strength 
imposed.    A  constriction  is  suggested  where  the  location  and  magnitude  of  these  minima 
appear  affected  by  the  initial  field  strength  as  shown  in  Figure  (8). 
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Figure  8:  Ion  Density  vs  Distance  -  E(0)  >  E^^ 
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Figure  (8)  suggests  the  extent  of  the  sheath  as  order  10("6),  in  agreement  with  the 
theory  which  requires  the  shielding  to  be  of  Debye  length:  AD=[e0kT0/q2no]H   But  only 
sufficiently  strong  initial  electric  field  intensities  are  able  to  push  the  ion  densities  past  their 
wells,  past  their  constriction  (Figure  (8)). 

In  comparison,  the  electrons  in  the  anode  sheath  behave  as  expected,  being  less 
influenced  by  electrostatic  repulsion  than  the  ions.  Figure  (9)  shows  that  for  very  low  initial 
field  strength,  the  electron  density  increases  monotonically,  but  only  very  slightly. 
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Figure  9:  Electron  Density  vs  Distance  -  E(0)  <,  E^,^ 

So,  it  is  expected  that  increasing  initial  field  strength  should  give  rise  to  increasing 
electron  density  -  this  is  the  general  case  (Figure  (10)).  The  electron  density  does  decrease 
very  slightly  before  the  location  of  ion  minimum  is  reached  as  illustrated  by  Figures  (10)-(1 1); 
Figure  (11)  represents  a  more  detailed  view  of  Figure  (10)  in  the  region  before  ion  minimum. 
Similar  to  the  ion  density,  only  sufficiently  large  initial  field  magnitudes  push  the  electron 
density  past  the  location  of  ion  well  (Figure  (10)).  From  a  physical  standpoint,  it  should  be 
remembered  that  both  species  are  affected  by  electrostatic  forces:  the  ions,  however,  are 
influenced  to  a  much  greater  extent  due  to  the  electrostatic  repulsion  required  when  like 
charges  are  sufficiently  close;  qualitatively,  these  are  the  observations  found  in  the  preceding 
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analysis  (Chapter  V). 
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Figure  10:  Electron  Density  vs  Distance  -  Near  and  Past  Density  Minimum  E(0)  >  E^b^m 
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Figure  1 1 :  Electron  Density  vs  Distance  -  E(0)  <  E^g^,,,  E(0)  >  Ethreshold  (Expanded  View) 
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If  the  initial  field  is  sufficiently  strong  to  push  the  ion  density  past  its  minimum,  the  minimum 
occurs  when  the  electric  field  curvature  changes:  that  is,  the  field  curvature  as  written  in  terms 

d2  I  ~\ 

of  its  potential     ^^\V)     changes  signs.     Shown  in  Figure  (12)  are  the  results  for 

dy2 


£(0)  =  5.08       with       ».(0)  =  108 


m 


,    «(0)=1012 


m 


Densities  &  Electric  Field  vs  Distance 
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Figure  12:  Densities  &  Electric  Field  vs  Distance  -  E(0)  >  EltacghoU=5.08 

Although  a  quasi-neutral  region  forms  when  the  initial  electric  field  is  sufficiently  strong  so 
as  to  push  the  ions  past  their  well,  the  number  density  near  the  plasma  boundary  is  not 

reached:  rc*  =0.15695       as     compared     to     the     maximum     density     achieved 

when  nt~n     where     n     is  on  the  order  of  10"3.  Since  nt~ne    requires     E    =  constant, 
something  else  is  needed  to  drive  the  system  toward  the  plasma  stagnation  -  some  form  of 
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transition.  To  summarize  the  effect  of  initial  conditions: 

It  appears  that  the  anode  fall  is  comprised  of  diverse  behavior,  in  part  determined  by 

the  initial  conditions  chosen:  there  appear  ranges  of  A«  =  \n.  -  »J     vs     E(0)     for  which 

both  equations  (7-9)  and  equations  (10)  naturally  form  a  region  of  quasi-neutrality  whereas 
for  other  ranges  of    Ah  =  \n.  -  nj      vs    £"(0)    ,  there  seems  insufficient  initial  energy  to 

push  the  ion  density  past  its  minimum.  To  gauge  this  effect,  equations  (10)  were  integrated 
again  but  now  E(0)  was  varied  from  0.01  to  20(1 04)  [V/m]  for  fixed  An  ranging  between  109 
to  1016  [1/m3].  Figure  (13)  shows  the  results  in  semi-logarithmic  format:  the  data  points  are 
plotted  for  which  E(0)  is  sufficient  to  drive  the  ion  density  past  the  minimum,  toward  the 
electron  density.  From  this  plot,  a  type  of  "cut-off'  region  seems  to  form:  to  the  left,  the  ions 
remain  locked  in  their  density  well,  to  the  right,  the  ions  are  driven  past  this  minimum  toward 
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Figure  13:  How  to  Begin  Selecting  Initial  Conditions  E  vs  An:  An=|ni-ne 
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For  the  region  of  initial  conditions  for  which     n,  ~n      is  reached,  an  interesting 

result  appears  when  the  ion  density  is  plotted  against  the  corresponding  ion  density:  a  type 
of  phase  portrait  for  the  ions.  Specifically,  a  multi-dimensional  homoclinic-behavior  results 
(see  Appendix  C),  a  characteristic  commonly  associated  with  energy  conservation:  for 
sufficient  initial  field,  energy  requires  a  balance  between  gains  and  losses,  between  ionization 
and  recombination.  So,  if  the  initial  conditions  are  such  that  the  ion  density  is  pushed  past  its 

well  toward  quasi-neutrality,  the  event  of  nt  ~  ne  closes  the  homoclinic  orbit  and  energy 
is  balanced:  a  comparison  between  the  case     E(0)  =  5.00    and     E(0)  =  4. 17    is  shown  in 


Figure  (14). 
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Figure  14:  Anode  Fall  with  Sufficient  Initial  Field  -  Phase  Portrait  for  Ions 


The  integrated  results  of  equations  (10)  not  only  provide  the  species  density  behavior, 
but  also  their  gradients.  As  such,  a  natural  question  arises  as  to  the  behavior  of  species  drift 
and  random  currents:  the  behavior  of  the  species  under  the  influences  of  an  electrostatic  field 
(potential)  as  compared  to  the  behavior  in  an  environment  free  of  potential. 
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2.  Drift  and  Random  Currents 

The  species  densities  are  not  constant  throughout  a  region  which  is  influenced  by  an 
electric  field;  a  uniform  electric  field  applied  to  a  cloud  of  constant  charge  q  will  produce  a 
directed  motion  along  with  a  variation  of  the  concentration.  From  this  principle,  the  species 
drift  currents  are  found  [Ref.  57:  pp.  53-55]: 


i.o  i.e 


Jo;  J  drift; .  ~*  Clnoieni,  e  Vavg  drift 


=  ±("»0      "/, 


M^f-^raKA.) 


0;."/,e  "-£ 


c^c- = -  n„   n,  A 


kT„ 


°i.»    i-e      E 


Equation  (13) 


where  the  Einstein  relation  for  mobility  has  been  applied.  To  obtain  the  random  current,  the 
thermal  velocity  is  used  such  that 


'o;  J  random; 


=  ± 


J )  ?^V««^-*(  4  j  Wo^e^ 


2kT„ 


m 


Equation  (14) 


'.e 


The  effect  of  varying  initial  conditions,  An  vs  E(0),  are  again  determined  in  similar  procedures 
as  for  the  species  densities  above  as  shown  next.  First,  the  ion  currents  are  considered: 

For  field  strengths      E(0 )<Ethreshold  =  0.04        ,  the  random  ion  current  increases, 

following  the  tendencies  of  the  ion  density  (see  Figure  (5)).  When  the  initial  field  is  increased 
such  that      E  ( 0 )  >  Ethreshold  =  0 .  04       ,  the  random  and  drift  currents  intersect  as  a  function 


of  the  rate  of  decrease  or  strength  of  initial  field.  Again,  the  random  current  (equation  (14)) 
follows  the  density  distribution  for  isothermal  conditions.  Figure  (15)  summarizes  these 
results  for  locations  "close"  to  the  electrode  surface. 
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Ion  Current  Density:  Random  &  Drift 

Varying  Initial  Electric  Field 
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Figure  15:  Random  &  Drift  Ion  Current  Density  vs  Distance  -  E(0)  <,  E^^,,, 

If  the  initial  electric  field  is  sufficiently  large  so  as  to  push  the  ions  past  their  density 
well,  the  drift  current  follows  accordingly  Furthermore,  the  random  current  follows  the  ion 
density  as  shown  in  Figure  (16). 
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Figure  16:  Ions  -  Drift  &  Random  Currents,  Density  &  Gradient  vs  Distance  -  E(0)=5.08 
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Of  note  is  illation  bolh_the  drift  and  jandom-cur-rents-seem  to  approach  a4imrtmg 
value- ai  the  realm  of  quasi-neutrality  is  approached;  perhaps  the_  saturation  current  as 
theorized  by  Child  and  Langmuir  [Ref.  58:  p.  238].  Next,  the  electron  currents  are 
considered: 

For  sufficiently  large  initial  field,  the  rapid  variation  in  the  electron  density  gradient 
(see  Figure  (1 1))  is  markedly  observed  in  the  corresponding  electron  drift  current  -  Figure 
(17).  Additionally,  the  electron  random  current  does  not  intersect  the  electron  drift  current, 
in  comparison  to  that  of  the  ions.  Instead,  the  random  electron  current  follows  the  tendencies 
of  the  electron  density  (Figure  (18))  where  the  negative  charge  of  the  electron  reverses  the 
magnitude  of  the  electron  density  plot  (see  Figure  (11)),  in  comparison  to  the  ion  results. 
Furthermore,  when  a  sufficiently  large  initial  field  is  applied,  the  electron  drift  current 
increases  past  the  point  (spatial  coordinate)  of  the  density  minimum,  diverging  spatially  from 
the  random  current  (Figure  (19)). 
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Figure  19:  Random  &  Drift  Electron  Current  Densities  vs  Distance  -  Near  and  After  Ion 
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B.         NUMERICALLY  MATCHING 

Although  a  quasi-neutral  region  forms  when  the  initial  electric  field  is  sufficiently 
strong  so  as  to  push  the  ions  past  their  well,  the  number  density  near  the  plasma  boundary  is 

not    reached:         n*  =  0.1 5695      as    compared    to    the    maximum    density    achieved 

when  ni~ne    where    n    is  on  the  order  of  10"3.  Since    nt~ne    requires     E    =  constant, 

something  else  is  needed  to  drive  the  system  toward  the  plasma  stagnation  -  some  form  of 
transition,  for  without  it,  the  species  densities  diverge  (Figure  (20,  21)). 
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Figure  20:  Anode  Fall  -  Without  Numeric  Matching 
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Figure  21:  Anode  Fall  -  Without  Numeric  Matching  (Expanded  View) 

In  steady-state  with  isothermal  conditions,  with  the  electric  field  constant  when 
ni  ~  ne     ,  there  remain  two  primary  phenomena  which  could  drive  the  density  toward  the 

plasma  equilibrium:  density  gradient  or  ambipolar  diffusion.  To  obtain  the  overall  effect,  first 
the  sheath  equations  as  given  by  set  (10)  are  matched  numerically  to  the  final  ambipolar 

equations,  (12)  using  for  the  matching  conditions      ni  ~  ne     that  result  from  equations  (10) 

provided  the  initial  field  is  sufficiently  strong  to  allow  the  ions  to  proceed  past  their  density 
well;  the  results  for  E(0)=6.1(104)  [V/m]  and  rii(0)=108  [1/m3]  with  ne(0)=1012  [1/m3]  are 
shown  in  Figure  (22). 
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Figure  22:  Anode  Fall  -  Numeric  Matching  without  Transition 
The  density  remains  below  the  plasma  boundary  stagnation  (Figures  (22,  23))- 
ambipolar  diffusion  is  insufficient  to  drive  the  density  toward  the  plasma  boundary.    The 
cosinusoidal  behavior  shown  is  determined  by  the  matching  conditions  imposed  by  the  sheath 
equations  (10)  upon  the  ambipolar  equation  (12)  -  Chapter  V. 
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Figure  23:  Anode  Fall  -  Numeric  Matching  without  Transition  (Expanded  View) 
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If,  on  the  other  hand,  a  full  transition  region  is  considered  in  which  the  density 
gradient  is  the  primary  source  to  drive  the  system,  then  the  plasma  boundary  density  is 
approached.  The  transition  equations  do  not,  however,  reflect  the  impact  of  the  density 
growth  on  free  diffusion  -  the  formation  of  ambipolar  diffusion.  Hence,  using  the  plasma 

boundary  stagnation,     n  *  =  0. 15695     ,  along  with  the  electric  field  intensity  at  the  onset  of 
ni~ne     as  initial  conditions  (or  matching  conditions)  to  equations  (12),  the  transition 


equations  (1 1)  are  numerically  matched  to  the  ambipolar  formulation,  set  (12).  Figure  (24, 
25)  illustrates  the  results  for  E(0)=6.1(104)  [V/m]  andni(0)=10?  [1/rrf  ]  withne(0)=1012  [1/irf]. 
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Figure  24:  Anode  Fall  -  Numeric  Matching  with  Transition 
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Figure  25:  Anode  Fall  -  Numeric  Matching  with  Transition  (Expanded  View) 

The  ambipolar  region  shown  in  Figures  (24,  25)  clearly  shows  the  very  sensitive 
nature  of  this  zone,  the  initial  conditions  to  the  transition  region  the  consequent  matching  to 
the  ambipolar  zone  determine  the  characteristics  of  this  latter  region  (see  Chapter  V,  Figure 
(4)).  Specifically,  the  behavior  of  the  anode  fall  region,  as  modeled  by  equations  (7-9)  and 
equations  (10-12),  is  governed  by  the  initial  conditions  imposed  to  the  transition  by  the 
sheath.  The  ultimate  condition  that  determines  the  anode  fall  characteristics,  however,  is  the 
initial  field  strength  imposed  upon  the  region:  if  the  field  is  sufficiently  strong,  the  ion  density 
is  pushed  past  its  well  toward  quasi-neutrality  -  these  matching  conditions,  for  the  transition 
equations,  determine  the  rate  of  density  increase  and  consequently  the  initial  conditions  for 
the  ambipolar  equations.  These  criteria  will  then  impose  the  behavior  for  the  ambipolar 
region,  whether  the  manifolds  of  the  plasma  boundary  stagnation  restrict  the  system 
characteristics,  or  whether  the  realm  of  instability  is  entered  (see  Chapter  V  and  Figure  (4)). 
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VH.  SUMMARY,  CONCLUSIONS  AND  RECOMMENDATIONS 

A.         SUMMARY 

Voltage  losses  at  electrode  boundaries,  surface  erosion  and  sheath  effects  play  an 
important  role  in  plasma  devices  and  must  be  controlled  in  designs  of  practical  interest. 
Particularly,  thermal  arcjets  and  MPD  accelerators  deposit  between  15%  and  80%  of  the 
input  power  to  the  anode.  This  presents  not  only  a  severe  performance  penalty,  but  also 
complicates  the  thermal  design  problems  since  the  heat  thus  generated  must  be  radiated  away 
from  the  thruster  surfaces.  [Ref.  7,  9]  A  natural  question  arises:  what  is  the  behavior  found 
in  the  anode  region?  Perhaps  a  restatement  of  the  qualitative  picture  of  the  anode,  as  offered 
by  Ingold  [Ref.  43:  pp. 59-63],  serves  as  an  answer: 

1.  The  anode  fall  voltage  is  on  the  order  of  the  ionization  potential  of  the  gas  since 
electrons  must  be  accelerated  to  an  energy  sufficiently  high  so  that  an  ion  current  for 
the  positive  column  is  provided. 

2.  The  anode  fall  region  is  determined  primarily  by  space-charge  and  secondly  by  the 
ionization  requirements. 

Moreover,  in  a  comprehensive  review  of  anode-spot  phenomena,  Miller  [Ref.  44] 
suggests  that  a  sudden  increase  in  the  voltage  noise  frequently  indicates  that  an  anode  foot- 
print could  be  forming  or  that  a  transition  into  the  anode-spot  mode  could  be  occurring. 
Further,  this  abrupt  change  in  voltage  has  been  associated  with  a  sudden  change  of  ion  density 
in  the  region  near  the  anode,  a  density  decrease  that  would  leave  the  local  electron  space 
charge  uncompensated  and  thus  produce  a  low-conductivity  region.  The  sudden  increase  ;n 
voltage  would  then  reflect  this  local  shortage  of  ions  -  an  ion  depletion  region  or  ion 
starvation  region.  Even-though  the  most  suitable  model  for  the  anode  region  is  still  actively 
debated,  the  general  idea  of  the  appearance  of  an  ion  deficiency  region  near  the  anode  as 
triggering  the  transition  into  the  anode-spot  mode  has  been  fairly  supported. 

So  to  expand  these  qualitative  portraits  of  the  anode  fall,  this  work  investigated  the 
nature  of  the  voltage  drops  in  the  vicinity  of  a  non-emitting,  positive  electrode.  The  selected 
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approach  involves  non-linear  analysis  techniques  of  the  continuum  governing  equations  for 

steady-state,  isothermal  conditions  in  one  dimension,  where  both  ionization  and  two-body 

recombination  processes  are  considered.   The  following  conclusions  and  observations  are 

offered: 

B.         CONCLUSIONS 

Ions,  created  through  electron  bombardment  of  neutral  atoms  throughout  the  anode 
fall  region,  are  repelled  electrostatically  from  the  positive  electrode  toward  a  location  of 
quasi-neutrality  at  the  sheath  boundary,  as  well  as  toward  an  unstable  stagnation  region  near 
the  plasma  interface.  These  equilibria  form  as  a  result  of  the  chemistry  present.  That  is, 
continuity  (equations  (9))  seems  to  require  flow  equilibrium  when  quasi-neutral  conditions 
are  initially  met  within  the  sheath  and  when  the  plasma  boundary  is  reached.  It  is 
recombination  that  establishes  the  latter  while  the  presence  of  ionization  stipulates  the  former. 
Fundamentally,  the  sheath,  as  formulated,  is  inherently  unstable  -  ions  are  driven  out  of  the 
system,  toward  the  negative  electrode.  The  stagnation  regions  were  obtained  analytically, 
through  non-linear  analysis  techniques;  to  test  these  results,  nitrogen  data  were  used  in  a 
numerical  algorithm  from  which  the  following  observations  are  made: 

At  the  electrode  wall,  in  its  immediate  vicinity,  at  least  two  major  processes  repel  the 
positive  ions  away  from  the  wall  in  the  isothermal,  steady-state  case:  the  initial  electric  field 
strength  and  the  initial  density  gradient.  If  the  field  strength  is  low  enough,  the  gradient 
dominates  the  system  and  the  ion  density  behaves  as  expected,  increasing  monotonically  for 
a  short  distance.  If  the  field  strength  is  high  enough,  the  gradient  influence  becomes  secondary 
in  nature,  the  field  then  driving  the  species  densities.  However,  the  low  initial  field  strength 
and  ion  density  gradient  are  insufficient  to  push  the  species  toward  quasi-neutrality,  when 
i\~ne.  On  the  other  hand,  when  the  initial  field  becomes  sufficiently  strong  as  compared  to 
the  species  density  gradient,  the  ions  are  pushed  toward  what  appears  as  a  density  trough. 

Furthermore,  for  certain  initial  values  of    An=\nj-ne\     vs  E{ 0 )   ,  the  ions  are  pushed  past 

their  well,  eventually  approaching  the  electron  density  in  magnitude  but  not  merging:  for 
sufficiently  strong  initial  field,  quasi-neutrality  forms  in  a  self-consistent  way.    Further,  if 
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conditions  permit  the  ions  to  traverse  their  density  well,  their  constriction,  then  at  the  moment 
of  quasi-neutrality,  when     w. ~ne    ,  there  appears  a  conservative  phenomenon  for  the  ions, 

where  energy  growth  is  matched  by  energy  loss  -  a  chemical  balance. 

At  the  onset  of  quasi-neutrality,  however,  the  ensuing  species  densities  do  not  match 
the  plasma  neutrality  requirements  (Sana).  Thus,  a  transition  region  must  drive  the  system 
toward  the  plasma  boundary:  for  steady-state  conditions,  the  primary  mechanism  is  the 
density  gradient.  When  the  species  density  becomes  large  enough,  mutual  Coulomb  fields 
affect  free  diffusion  and  ambipolar  diffusion  results.  The  ambipolar  field  then  continues  to 
drive  the  ions  away  from  the  positive  plate;  the  initial  conditions  to  the  transition  region,  the 
consequent  matching  to  the  ambipolar  zone,  determine  the  characteristics  of  this  latter  region. 
Specifically,  the  behavior  of  the  anode  fall  region,  as  modeled  by  equations  (7-9)  and 
equations  (10-12),  is  governed  by  the  initial  conditions  imposed  to  the  transition  by  the 
sheath.  The  ultimate  condition  that  determines  the  anode  fall  characteristics,  however,  is  the 
initial  field  strength  imposed  upon  the  region:  if  the  field  is  sufficiently  strong,  the  ion  density 
is  pushed  past  its  well  toward  quasi-neutrality  -  these  matching  conditions  for  the  transition 
equations  determine  the  rate  of  density  increase  and  consequently  the  initial  conditions  for  the 
ambipolar  equations.  These  criteria  then  determine  whether  the  manifolds  of  the  plasma 
boundary  stagnation  restrict  the  system,  or  whether  the  realm  of  instability  is  entered.  As 
such,  any  numeric  solutions  will  most  likely  be  challenging,  a  result  observed  when,  using 
nitrogen  data,  the  various  anode  fall  regions  are  numerically  matched. 

Both  the  analysis  and  the  numerical  results  suggest  that  the  anode  fall  model  used  in 
this  work  does  not  represent  sufficient  dissipation,  preventing  the  system  from  adjusting  itself, 
to  adhere,  to  the  "stable  manifold"  of  the  stagnation  region  near  the  plasma  boundary.  The 
model  is  most  likely  incomplete  in  its  physical  representation:  either  temperature  gradients, 
diffusion  written  in  terms  of  the  electric  field,  full  time  dependencies  or  quasi-steady 
formulation,  magnetic  field  effects  and/or  three-body  recombination  effects  may  provide  the 
necessary  dissipation. 


67 


C.         RECOMMENDATIONS 

To  improve  upon  the  results  obtained  in  this  work,  further  considerations  should  be 
given  to: 

1 .  Species  temperature  variation  throughout  the  fall  region,  specifically  in  the  sheath. 

2.  Diffusion  perhaps  written  as  a  function  of  the  electric  field  instead  of  as  assumed 
constant  for  each  species. 

3.  Full  time-dependent  and/or  quasi-steady  formulations  for  the  anode  fall. 

4.  The  effect  of  three-body  recombination  on  the  stagnation  regions  found,  as  well 
as  the  influences  of  magnetic  fields  on  the  system  behavior. 

This  work  presents  a  broad  overview  of  many  subjects,  many  of  which  really  need 
further  investigation,  among  them: 

1 .  Analytically  matching  the  second  order  sheath  equations  to  the  transition  set  and 
to  the  ambipolar  equation;  essentially,  the  analytic  analog  of  the  numeric  matching 
ideas  presented  in  this  work.  The  results  of  such  mathematical  treatment  might  reveal 
that,  in  fact,  more  physics  is  needed  in  equations  (7-9). 

2.  Further  investigation  of  the  sensitivity  to  initial  conditions  for  the  ion  density  well 
as  well  as  quasi-neutral  zone  formation;  a  wider  range  of  initial  conditions  should  be 
tested  in  the  hopes  of  finding  some  better  relation  between  An  versus  E  which  would 
be  of  aid  in  any  numeric  endeavor.  Such  analysis  should  also  focus  on  a  variety  of 
gases,  not  be  limited  to  nitrogen;  possible  candidates  include  xenon  and  argon. 

3.  Generation  of  a  current  versus  voltage  curve,  I-V  curve,  by  integrating  the  electric 
field  data  and  using  the  current  calculations  given  in  this  work.  Does  the  result  match 
the  literature? 
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APPENDIX  A.  NON-DIMENSIONALIZATION 

Mathematical  derivations  are  frequently  burdened  by  complicated  functions  of  the 
constants  in  the  problems.  Before  the  mathematical  analysis  is  carried  out,  a  preliminary 
dimensional  analysis  may  not  only  reveal  the  ways  in  which  some  of  the  constants  enter  into 
the  final  solution,  but  also  unveil  significant  dimensionless  products  of  a  problem.  Then  the 
original  differential  equations  may  be  expressed  in  terms  of  dimensionless  notation,  a  method 
particularly  useful  since  model  laws  are  usually  derived  from  the  differential  equations  that 
govern  phenomena.  [Ref.  59:  p.  144]  Dimensional  analysis  is  a  process  in  as  much  influenced 
by  art  as  it  is  by  scientific  methods:  two  principles  were  applied  to  this  work  -  the 
Buckingham  Pi  Theorem  and  Fractional  Analysis  coupled  with  the  physical  requirements 
brought  by  the  governing  equations  themselves. 
A.         BUCKINGHAM  PI  THEOREM 

The  Pi  Theorem  may  be  stated  as  "The  number  of  dimensionless  products  in  a 
complete  set  is  equal  to  the  total  number  of  variables  minus  the  rank  of  their  dimensional 
matrix."  [Ref.  59:  p.  31]  In  using  the  Pi  Theorem,  the  following  conditions  should  be  fulfilled: 
[Ref  60:  p.  20] 

1.  The  list  of  dimensional  parameters  must  contain  all  of  the  parameters  of  physical 
significance  including  all  independent  parameters  and  one  dependent  parameter. 

2.  The  non-dimensional  pi's  as  finally  composed  should  contain,  at  least  once,  each 
of  the  parameters  in  the  original  list. 

3.  The  list  of  dimensions  used  to  compose  the  physical  parameters  must  be 
independent,  or  else  provision  must  be  made  to  compensate  for  the  redundancy. 

There  are  infinitely  many  different  complete  sets  of  dimensionless  products  that  can 
be  formed  from  a  given  set  of  variables.  Insofar  as  Buckingham's  theorem  is  concerned,  any 
such  complete  set  is  admissible,  with  some  sets  more  useful  than  others.  So  how  may  a 
complete  set  of  dimensionless  products  be  most  advantageously  selected  at  the  outset? 
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1.         Arrangement  of  Variables  and  Terminology 

As  Buckingham  has  pointed  out,  the  maximum  amount  of  experimental  control  over 
the  dimensionless  variables  is  obtained  if  the  dependent  variable  does  not  occur  in  more  than 
one  dimensionless  product,  a  product  referred  to  as  the  "dependent  dimensionless  variable 
or  pi"  [Ref.  59].  Using  a  dimensional  matrix  as  outlined  in  reference  59,  the  preceding 
condition  will  be  realized,  as  nearly  as  possible,  if  in  the  dimensional  matrix  the  first  variable 
is  the  dependent  variable,  the  second  variable  is  that  which  is  easiest  to  regulate 
experimentally,  the  third  variable  is  that  which  is  next  easiest  to  regulate  experimentally,  and 
so  on.  [Ref.  59:  p.  39]  Consequently,  the  variables  describing  equations  (7-9)  are  arranged 
as  follows:  E=xx,  n;=x2,  n^,  je=x4,  y=x5,  De=X6,  v=x7,  j;=xg,  D^,  k=x10,  T0=xn,  cc2=x12, 
e0=x13,  and  q=xu  where  E  represents  the  electric  field,  r^e  characterize  the  ion  and  electron 
number  densities,  respectively,  j^e  are  the  ion  and  electron  current  densities,  respectively,  y  is 
the  geometry  coordinate  length,  D^e  characterize  the  ion  and  electron  difrusivities, 
respectively,  and  v  is  the  ionization  coefficient.  Universal  constants  and  those  variables  held 
constant  are  entered  last  in  the  matrix  since  these  are  not  deemed  "experimentally 
controllable":  k  is  the  Boltzmann  constant,  a2  the  two-body  recombination  coefficient,  and 
e0  the  permittivity  of  free  space  with  q  as  the  charge.  The  entries  of  the  matrix  are  next  made: 
the  dependent  variable,  i.e.,  the  first  entry  in  the  matrix,  is  the  electric  field,  since  in  terms  of 
voltage  applied,  it  governs  the  space  charge  development.  The  second  and  third  entries  into 
this  matrix  represent  the  number  density  distributions  and  govern  the  current  development. 
All  other  variable  entries  are  made  arbitrarily,  the  only  requirement  being  that  arrangement 
which  gives  a  linearly  independent  solution  for  x10. .  .x14.  To  complete  the  dimensional  matrix, 
each  variable  and  constant  is  written  in  terms  of  the  fundamental  units  of  mass  (M),  length 
(L),  time  (t),  charge  (Q)  and  temperature  (0);  the  matrix  outlining  the  coefficients  of  these 
units  for  each  variable  becomes: 
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Table  7:  Variables  in  Terms  of  Fundamental  Units 
There  will  be  14-5=9  dimensionless  pi's  (number  of  variables  minus  the  number  of  units 
necessary  to  describe  them).  In  fact  the  rank  of  this  matrix  is  five  which  not  only  reflects  the 
number  of  fundamental  units  necessary,  but  also  the  five  constants  of  the  problem:  k,  T0,  q, 
ew  (%.  Solving  for  the  variables  xj  for  i=10..  14,  in  terms  of  the  variables  Xj  for  j=l  .9  results 
in  the  homogenous,  linear,  algebraic  equations.  By  taking  this  equation  set  into  matrix  format 
and  solving  for  each  variable  X;  where  i=1..9,  the  homogeneous,  linearly  independent  non- 
dimensional  pi's  can  be  determined: 


MO 
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Table  8:  Buckingham  Pi's  -  Matrix  Formulation 
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Non-dimensional 
Variables 

Relationship  Between  Non-dimensional 
Variables,  Buckingham  Pi's  and 
Undisturbed  Plasma  Values 
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Table  9:  Relationship  Between  Non-dimensional  Variables,  Buckingham  Pi's  and 

Undisturbed  Plasma  Values 


When  these  pi's  are  inserted  into  the  appropriate  governing  equations,  each  term  in  the 
equation  will  have  a  non-dimensional  parameter  affecting  it.  For  example,  the  dimensional 

Gauss       —  =—(»*- n.)      written  with  Buckingham  Pi ' s  gives 
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It  is  the  constants  and  pi's  within  the  brackets  that  form  the  non-dimensional  parameters. 
2.         Pi  Theorem  -  Results  and  Conclusions 

Although  the  Pi  Theorem  forms  a  good  framework  for  discussing  the  nature  of  units, 
dimensions  and  related  topics,  the  theorem,  when  used  as  the  only  means  toward  analysis, 
suffers  from  several  deficiencies: 

1 .  No  direct  means  for  finding  the  pertinent  pi's  is  available  and  dimensional  analysis 
itself  provides  little  framework  for  incorporating  relevant  physical  information. 

2.  The  theorem  alone  does  not  provide  conditions  under  which  one  or  more  pi's  can 
be  neglected  or  for  which  sets  of  dimensionless  groups  may  be  combined. 

So  another  method  must  be  used  to  make  the  problem  non-dimensional. 
B.         FRACTIONAL  ANALYSIS 

There  is  another  common  method  of  dimensional  analysis  that  is  intrinsically  the  same 
as  the  foregoing  technique  but  differs  from  it  in  that  the  differential  equations  are  expressed 
in  dimensionless  forms  by  applying,  among  others,  a  characteristic  length  L  and/or  a 

—    x      ~    t 
characteristic  time  period  x,  such  that     x  =  — ,    T=—  ;  that  is,  all  variables  are  non- 

L  x 

dimensionalized  by  introducing  appropriate  characteristic  values  -  fractional  analysis.  [Ref 
59]  A  fractional  analysis  is  any  procedure  for  obtaining  some  information  about  the  answer 
to  a  problem  in  the  absence  of  methods  or  time  for  rinding  a  complete  solution,  a  method 
from  which  an  approximate  analytical  or  numerical  solution  can  be  obtained  [Ref.  60]. 
Specifically: 
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1.  What  is  the  physical  meaning  of  each  of  the  governing  parameters  and  variable? 
What  are  the  qualitative  effects  of  an  increase  or  decrease  in  any  given  parameter  or 
variable? 

2.  Can  we  find  the  conditions  under  which  the  effects  of  certain  parameters  can  be 
neglected  either  in  a  given  region  or  for  a  particular  problem?  If  so,  does  this  lead  to 
governing  equations  that  are  more  tractable  so  that  they  can  be  solved  even  though 
the  general  equations  cannot  be  solved? 

3.  Are  there  any  combinations  of  two  or  more  non-dimensional  parameters  or 
variables,  or  transformations  of  variables,  which  lead  to  fewer  independent  quantities 
or  which  simplify  the  correlations  achieved? 

To  answer  these  questions,  the  confines  of  traditional  dimensional  analysis,  as  a  tool 
by  itself  are  exceeded  and  other  methods  must  be  employed.  One  such  method  incorporates 
dimensional  analysis  with  governing  equations  as  advocated  by  Sedov  [Ref.  61].  Sedov  has 
clearly  shown  that  important  information  can  be  obtained  by  simultaneous  use  of  the 
governing  equations  together  with  dimensional  analysis.  Consequently  a  systematic 
methodology  is  needed  to  obtain  information  directly  from  the  governing  equations  without 
actually  solving  them,  a  methodology  involving  three  primary  ideas:  [Ref.  60:  p.  66-70] 

1 .  Normalization  based  on  governing  equations;  making  equations  and  conditions 
non-dimensional  in  terms  of  non-dimensional  variables  of  standard  magnitude. 

2.  Absorption  of  parameters. 

3.  Combination  of  variables. 

To  carry  out  such  procedures,  a  clear  understanding  of  the  physical  information 
inherent  in  the  equations  is  essential.  Of  equal  importance  is  using  a  standard  procedures  for 
transforming  the  variables  to  non-dimensional  form  and  standard  magnitude;  it  is  not 
sufficient  to  make  governing  equations  non-dimensional  in  any  arbitrary  way.  So  of  the  three 
processes  stated  above,  normalization  is  the  most  important.  [Ref.  60] 

Here,  normalization  is  defined  as  making  the  governing  equations  and  conditions  non- 
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dimensional  in  terms  of  non-dimensional  variables  of  standard  magnitude.  To  carry  out  a 
normalization,  two  steps  are  usually  required:  1)  making  all  the  variables  non-dimensional  in 
terms  of  the  appropriate  scales  of  the  problem;  2)  dividing  through  the  equation  by  the 
coefficient  of  one  term  to  make  the  equation  dimensionless  (unit-free)  term  by  term.  Further, 
the  method  for  choosing  the  scaling  is  critical  as  this  choice  will  permutate  throughout  the 
ensuing  analysis.  To  this  end,  the  following  procedure  is  used  (after  Kline): 

1 .  Define  all  dependent  non-dimensional  variables  so  that  they  are  approximately 
unity  over  a  finite  distance  and  nowhere  exceed  approximately  unity  in  the  domain  of 
concern;  in  this  work,  the  undisturbed  plasma  values  are  used,  denoted  by  a  "o" 
subscript. 

2.  Define  all  independent  non-dimensional  variables  so  that  their  increment  is 
approximately  unity  over  the  same  domain  of  concern  (0  to  1,  1  to  2,  etc,  in  the  new 
variables). 

As  a  result,  the  dimensionless  groups  formed  are  usually  composed  from  the  boundary 
conditions;  from  the  characterizing  sizes  or  scales  of  the  body;  and  from  the  physical 
constraints  of  the  original  equation  such  as  system  properties,  physical  constants  or  both. 
In  general,  it  is  not  known  in  advance,  which  choice  of  pi's  gives  the  most  useful  set  of 
parameters  but  it  is  obvious  that  the  fewer  pi's  required  to  specify  the  behavior,  the  more 
useful  the  result  will  be.  So,  how  to  reduce  the  number  of  pi's  and  consequently  parameters? 
[Ref  60:  pp.  75-94]  The  choice  of  length  or  time-scales  coupled  with  physical  insight  are  the 
first  step  in  this  effort,  an  effort  leading  to  a  "modified  fractional  analysis." 
C.         MODIFIED  FRACTIONAL  ANALYSIS:  ANODE  FALL 

1.         Non-dimensional  Parameters 

The  derivative  term  found  in  the  governing  equations  for  the  anode  fall,  equations 
(7-9),  is  important  since  without  which  the  conditions  at  the  plasma  boundary  could  not  be 
satisfied.  So  to  make  all  the  variables  non-dimensional  in  terms  of  the  appropriate  scale,  the 
equations  are  divided  through  by  the  coefficients  of  the  derivative  terms  to  form  the  necessary 
system  parameters  where  both  Buckingham  Pi's  and  undisturbed  plasma  values  were  used  to 
normalize  the  variables.  Following  is  a  list  of  the  results  obtained: 
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Governing  Equations 
[Equations  (7-9)] 


Non-dimensional 

Parameters  via 

Buckingham  Pi's 
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Parameters  via  Undisturbed 

Plasma  Values 
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Non-dimensional  Parameters 


2.         Estimate  of  the  Length  Scales 

In  the  above,  the  variables  where  normalized  using  undisturbed  plasma  values  since 
the  Buckingham's  Theorem  does  not  reflect  the  physics  at  hand:  for  example,  the  process  of 
ionization  affects  the  number  density  distribution  in  the  sheath  resulting  in  electric  field 
variations  and  consequently  influencing  the  current  distribution.  It  would  be  difficult  for  one 
length  scale  to  capture  each  of  these  behaviors.  Consequently  physical  reasoning  must  be  the 
foundation  from  which  to  estimate  the  length  scales  of  the  problem. 

The  neutral  plasma  is  shielded  from  the  space-charge  region  through  coulombic 
effects,  effects  which  should  drive  the  length  scales  affecting  the  associated  momentum 
exchange  processes.  The  first  term  on  the  right-hand  side  of  equations  (7a)  and  (7b)  model 
these  effects  and  are  used  to  estimate  the  length  scale  affecting  the  number  density 


distribution: 


L  ,  L 


Specifically, 
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~^        Equation  (7a) 
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Equation  (7b) 


where  [a]  and  [c]  are  the  coefficients  for  the  coulombic  terms  summarized  in  Table  2  (main 
body  of  this  work).  Through  the  process  of  non-dimensionalization,  the  parameters  are  O(l) 


kT  kT 

and  the  length  scales  become:       L  =  — -,   L  =  — - 


Similarly,  the  length  scale  for  the 


electric  field  comes  strictly  from  Gauss,  equation  (8),  where  either  term  can  be  used  to 
estimate  the  electric  field  length  scale  since  the  undisturbed  plasma  value  for  n0  is  the  same 
for  both  ions  and  electrons: 
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dy 
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Equation  (8) 


Using  the  non-dimensionalization  result  that  all  parameters  are  0(1),       LE 


4"0 


where 


the  parameters  [e]  and  [f]  are  again  given  in  Table  2. 

Within  the  sheath,  the  ionization  process  outweighs  the  recombination  loss,  so  in 
equations  (9a)  and  (9b)  the  ionization  term,  the  first  term  on  the  right-hand  side  of  the 
equations,  is  used  to  estimate  the  current  density  length  scales: 

d(L) 


dy_ 
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\ne  + 1  h 


nine      Equation  (9a) 
r  Equation  (9b) 


That  is,  the  non-dimensional  parameters  [g]  and  [i]  (Table  2)  estimate  the  length  scales  for 
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the  appropriate  current  density:         L= — — ,   L.= — - —        .  Essentially,  the  current  scale 


Jt    qvnn         h     q\nn 

*     o    o,  *■     o    o. 


reflects  the  influence  of  ionization. 

3.  Dimensionless  Products  -  Groupings 

Occasionally  transforming  dimensionless  products  to  achieve  greater  control  of  the 
variables  is  desired  or  may  form  a  useful  grouping:  the  various  results  formed  through  non- 
dimensionalizing  the  Navier-Stokes  equations  are  examples.  The  following  combination  of 
parameters  have  been  found  from  equations  (10-12): 
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Parameter  Groupings 
from  Equations  (10-12) 

Groupings  in  Terms  of  Problem  Constants 
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Table  10:  Non-dimensional  Parameter  Groupings 
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4.  Summary 

Fractional  analysis  coupled  with  the  use  of  the  governing  equations  indicates  that 
potentially  five  length  scales  are  present,  one  for  each  of  the  governing  equations;  however, 
by  using  undisturbed  plasma  values  for  the  normalization,  three  length  scales  result:  a  scale 
reflecting  the  conservation  of  momentum,  one  for  the  variation  of  the  electric  field,  and  the 
third  for  charge  conservation.  Of  these  three  lengths,  the  electric  field  scale  is  the  smallest, 
indicating  the  extent  of  the  sheath  as  the  space-charge  region  defines  this  region.  The  number 
density  scale,  on  the  other  hand,  acts  over  the  entire  range  of  the  anode  fall  and  so  should  be 
larger  in  magnitude  than  that  for  the  electric  field  scale.  The  largest  of  the  length  scales  is  the 
current  scale  since  the  current  follows  the  variation  in  species  density  according  to  js=qnsVs 
where  s=i,e: 
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Length  Scales  via  Fractional  Analysis 
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APPENDIX  B.  IONIZATION  AND  TWO-BODY  RECOMBINATION  DATA 
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Figure  26:  Ionization  Coefficient  v/N  as 
Function  of  E/N  for  N2  [after  Ref.  41] 
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Figure  27:  Two-Body  Recombination 
Coefficient  as  Function  of  E/N  for  N2 
[Ref  40] 


Note:  1  Td=10'21  [Vm2] 
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APPENDIX  C.  NON-LINEAR  DYNAMICS  -  AN  OVERVIEW 

When  exact  closed-form  solutions  of  a  physical  system  as  modeled  by  appropriate 
differential  equations  are  difficult  to  obtain,  or  when  the  exact  solution  is  too  complicated  to 
be  useful,  then  the  first  step  toward  an  approximate  solution  is  local  analysis.  The  results  of 
such  analysis  are  valid  only  in  a  sufficiently  small  neighborhood  of  a  point,  so  ultimately  a 
uniform  approximation  to  the  behavior  of  the  solution  over  an  entire  interval  may  be  found 
by  piecing  together  regions  where  the  local  behavior  is  known  through  global  analysis 
techniques  involving  perturbative  and  asymptotic  methods.  The  following  is  a  brief  overview 
and  adaptation  of  references  55  and  54  as  well  as  class  notes,  reference  62. 
A.         LOCAL  ANALYSIS:  THE  BASICS 

If  the  general  system  given  by 


is  visualized  as  trajectories  flowing  through  an  n-dimensional  phase  or  velocity  space  with 
coordinates  (x,,..,^),  differential  equations  can  be  interpreted  geometrically  as  vector  fields, 
where  now  the  velocity  of  the  flow,  the  rate  of  change  of  the  system  behavior,  can  be  plotted 
and  the  corresponding  equilibria  (stagnation  or  fixed  points)  of  the  system  determined  where 
the  flow  velocity  is  zero:  the  fixed  points  are  calculated  when  the  first  variation  tends  to 
zero  -  calculus  of  variations.  Linearization  techniques,  local  Taylor  series  expansion  about 
the  stagnation,  can  now  be  used  to  find  the  stability  about  these  points.  Such  analysis 
contains  qualitative  information  about  the  dynamical  behavior  of  the  system  leading  to  an 
overall  explanation  as  to  "why  things  happen  the  way  they  do." 
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In  general,  the  behavior  of  a  one-dimensional  dynamical  system        x  =f(x)        is 

considered  as  a  fluid  flowing  along  the  real  line  with  local  velocity/ (x).  This  imaginary  fluid 
is  referred  to  as  the  phase  fluid  and  the  real  line  is  the  phase  space.  In  higher  dimensions,  a 

similar  situation  occurs.  The  solution  to       x  =/(Xj,...,xw)    V    i=l,2,  •,«       given  an 

arbitrary  initial  condition  Xj,  is  then  represented  by  the  trajectory  of  the  flow  in  /^-dimensional 
space  as  determined  by  the  function  under  investigation.  The  collection  of  these  trajectories 
are  referred  to  as  a  phase  portrait  whose  appearance  is  controlled  by  the  fixed  points 
(stagnation,  equilibria  or  zeros)  x*  as  defined  by/(x*)=0.  This  fixed  point  is  defined  as  stable 
if  all  sufficiently  small  disturbances  away  from  it  dampen  out  as  t-°°.  Conversely,  unstable 
equilibria  are  those  for  which  disturbances  grow  in  time. 
B.         LINEAR  STABILITY  ANALYSIS 

An  n-th  order  differential  equation  can  be  expressed  as  a  system  of  n  first-order 
equations  by  appropriate  changes  in  the  variables.  To  determine  the  rate  of  decay  to  a  stable 
fixed  point,  linearization  techniques  about  the  equilibria  are  applied.  For  example,  let  x*  be 
a  fixed  point  and  let  r)(t)=x(t)-x*  be  a  small  perturbation  away  from  x*.  To  see  whether  the 
perturbation     grows     or     decays,     a     differential     equation     for     r\     is     formed: 

f|  =  — (x-x*)=x     since  x*  remains  constant  by  definition  of  stagnation  points. 
dt 


Consequently, 

i]  =  ^-(x-x*)=x=f(x)=f(x*+r])    where     f(x*  +  r\)»f(x*)  +  rtf'(x')  +  0(i\2) 
at 


84 


Noting  that/(x*)=0,      r\  =  r\f,(x*)+0(r\2)      .  Moreover,  if/  '(x*)*0,  then 

f|  =  r\f\x  *)        is  the  linearization  about  x*  where  the  perturbation  r|(t)  grows 

exponentially  if /'(x*)>0  and  decays  exponentially  provided/  '(x*)<0.  In  essence,  the  slope 
or  the  rate  of  change  of/'(x*)  at  the  fixed  point  determines  the  stability  of  the  stagnation. 
If  the  dimension  of  the  dynamical  system  is  more  than  one  dimensional,  the  procedure 
for  stability  analysis  is  similar  but  now  a  multi-variable  Taylor  series  expansion  is  used  for  the 
linearization  technique.  Consider  a  two-dimensional  non-linear  system  where  the  general 
form  of  a  vector  field  on  the  phase  plane  is  given  by 

X2~J2^Xl,X2^ 

with/  and/  as  given  functions.  In  vector  notation,  x=(x!,x2)  and/(x)  =  \fx  (x)/  (x)].  Here, 
x  represents  a  point  in  the  phase  plane  and    x    is  the  velocity  vector  at  that  point.   By 

flowing  along  the  vector  field,  a  phase  point  traces  out  a  solution  x(t),  corresponding  to  a 
trajectory  winding  through  the  phase  plane.  Essentially  the  entire  phase  plane  is  filled  with 
trajectories  since  each  point  can  play  the  role  of  an  initial  condition.  For  nonlinear  systems, 
however,  it  remains  a  challenging  task  to  find  these  trajectories,  the  solutions,  analytically. 
Even  when  explicit  solutions  are  available,  their  form  is  often  intractable  and  provide  limited 
insight:  the  importance  of  a  qualitative  analysis  to  understand  the  characteristics,  the 
dynamical  behavior,  of  a  system  becomes  apparent. 
Consider  the  system 

*=f(x,y) 

y=g(*,y) 
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and  suppose  that  (x*,y*)  represents  an  equilibrium  position,  i.e.,/(x*,y*)=0  and  g  (x*,y*)=0. 
To  begin  linearization,  let  w=x-x*  and  let  w=y-y*  be  components  of  small  perturbations  away 
from  the  stagnation.  To  see  whether  the  disturbance  grows  or  decays,  a  differential  equation 
for  u  and  w  must  be  found  similarly  to  the  one-dimensional  method  above  (noting  that  x*  and 
y*  remain  constant  by  definition  of  equilibrium): 


and 


u = x  =f(x  *  +uy  *  +w) 

dx  dy 

=  u—(f)+w  —  (f)  +  0{u2,w2,uw) 

ox  dy 


w  =y =g(x  *  +UJ  *  +w) 

=g(x*J*)+u  —  (g)+w  —  (g)  +  0(u2,w2,uw) 

dx  dy 

=  u  —  (g)+w  —  {g)  +  0{u2,w2,uw) 

dx  dy 


where  partial  derivatives  in  the  preceding  are  all  evaluated  at  (x*,y*).  That  is,  the  effects  of 
the  perturbations  (u,w)  evolve  according  to  the  linearized  system: 


dx  dy 

—(g)    —(g) 

dx  dy 


+  0(u2,w2,uw) 


(*V) 


The  matrix  of  partial  derivatives  evaluated  at  the  stagnation  (x*,y*)  is  termed  the  Jacobian 
matrix,  the  multi-variable  analogue  of  the  derivative/'(x*)  developed  for  the  one-dimensional 
system.  Furthermore,  if  the  quadratic  terms  found  in  0(u2,w2,uw)  are  sufficiently  small,  then 
the  behavior  of  the  linearized  system  can  be  analyzed  through  study  of  the  Jacobian  matrix 
alone. 
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C.         CLASSIFICATIONS  OF  FIXED  POINTS 

To  classify  all  possible  phase  portraits  that  can  occur  given  an  arbitrary  Jacobian 
matrix,  it  should  be  noted  that  the  coordinate  axes  play  a  critical  geometric  role;  they 
determine  the  direction  of  the  flow  as  t-=b»,  containing  straight-line  trajectories.  A  trajectory 
starting  on  one  of  these  axes  stays  on  that  axis  forever  and  exhibits  simple  exponential  growth 
or  decay  along  it.  For  the  general  case,  however,  the  analog  of  these  straight-line  trajectories 

has  to  be  found,  so  trajectories  of  the  form  x(t)  =  eXtV    ,  where  V>0  is  some  fixed  vector, 

are  sought  with  A  as  a  growth  rate,  also  to  be  determined.  If  such  solutions  exist,  they 
correspond  to  exponential  motion  along  the  line  spanned  by  the  vector  V.  The  conditions  on 
X  and  V  are  determined  by  the  eigenvalues  and  corresponding  eigenvectors  of  the  Jacobian 
matrix,  respectively.  That  is,  the  characteristics  of  the  linearized  system  are  determined  by 
the  nature  of  the  Jacobian  eigenvalues. 

1.  Robustness  of  Linearization 

The  stability  of  a  fixed  point  can  be  classified  in  terms  of  robust  or  marginal  cases. 
Robust  cases  are  those  for  which  the  stability  does  not  change  under  the  influence  of  small 
non-linear  terms  such  as  repellers  (sources),  attractors  (sinks)  or  saddles.  Specifically, 
repellers  are  those  fixed  points  for  which  the  Jacobian  displays  positive  real  parts  for  all 
eigenvalues,  attractors  are  those  characteristics  which  display  negative  real  parts  for  all 
eigenvalues.  Saddles  occur  when  one  eigenvalue  is  real  and  positive  while  the  other  is  real 
and  negative.  Marginal  stability  results  when  both  eigenvalues  are  imaginary  (centers)  or 
when  at  least  one  of  the  eigenvalues  is  zero.  For  marginal  stability,  then,  at  least  one  of 
eigenvalues  satisfies  Re(A)=0  whereas  for  robust  stability,  Re(A)*0;  that  is,  for  marginal 
stability,  the  trace  of  the  Jacobian  matrix  is  zero  or  the  flow  divergence  is  zero,  implying  that 
the  flow  volume  is  preserved. 

2.  Real,  Distinct  Eigenvalues 

In  general,  trajectories  approach  the  origin  of  the  phase  portrait  tangent  to  the  slow 
eigendirection,  the  direction  spanned  by  the  eigenvector  with  the  smaller  \k\.  If  |A2|<|Ai| ,  then 
the  trajectories  approach  A2  as  t-«>: 
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Fast 


y     \ 


W_ 


j  Slow 


Real,  Distinct  Eigenvalues  - 

0<|X2|<|Aj|  and  Xu  X2<0 

"Node" 

If  all  eigenvalues  are  positive  and  distinct,  then  the  behavior  locally  about  the  respective 
stagnation  is  termed  "unstable  node",  whereas  if  all  eigenvalues  are  negative  and  distinct,  then 
the  dynamics  about  the  fixed  point  is  called  "stable  node". 

What  happens  when  one  of  the  eigenvalues  is  positive  while  the  other  is  negative? 
If  XX),  the  corresponding  eigensolution  grows  exponentially;  if  X<0,  the  eigensolution  decays 
exponentially.  The  stable  manifold  corresponds  to  the  line  spanned  by  the  eigenvector  formed 
by  A<0,  whereas  the  unstable  manifold  corresponds  to  the  line  spanned  by  the  eigenvector 
formed  by  A>0.  The  combination  of  such  behavior  forms  the  "saddle"  where  trajectories  are 
attracted  toward  the  stable  manifold  and  repelled  as  the  unstable  manifold  is  approached: 


Unstable 


Real,  Distinct  Eigenvalues 
X,<0,  X2>0  "Saddle" 
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3.         Complex  Eigenvalues 

The  eigenvalues  of  the  Jacobian  determine  the  local  dynamical  behavior  of  the  system: 

the  eigensolution       x(t)=cle  ltV1+c2e^'  V2        is  affected  by  the  eigenvalues  according 

to  the  behavior  of  e^ .  If  Ai;2=a±ico,  then  the  system  is  influenced  by  e^e?"  where  a  is  the  real 
component  of  the  complex  eigenvalue  (damping)  and  go  is  the  imaginary  component,  o)*0. 
By  Euler's  formula,  e^*=cos(ot)+i[sin(cDt)]  and  so  x(t)  becomes  a  cosinusoidal  function  with 
exponential  growth,  instability,  according  to  eat  if  a>0,  or  with  exponential  decay,  stability, 
if  a<0.  Such  behavior  is  termed  "spiral",  unstable  or  stable  depending  on  whether  energy  is 
supplied  to  the  flow  or  taken  from  the  flow. 


Complex  Eigenvalues: 
A12=a±icot  "Spiral" 

I£  on  the  other  hand,  the  eigenvalues  are  purely  imaginary,  a=0,  then  all  solutions  are 
periodic  with  period  T=27t/cl)  where  the  oscillations  have  constant  amplitude. 


Imaginary  Eigenvalues: 
A12=±io)t  "Center" 
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Centers,  or  systems  with  purely  imaginary  eigenvalues,  occur  generally  where  energy 
is  conserved,  corresponding  to  a  state  of  neutrally  stable  equilibrium,  or  a  position  of 
minimum  energy,  although  the  linearization  represented  is  not  robust.  Small  non-linear  terms 
can  change  the  characteristics  from  periodicity  to  exponential  growth  or  decay  (see  below). 

4.  Real,  Equal  Eigenvalues 

Suppose  that  X1=X2=X,  the  eigenvalues  of  the  system  are  equal.  Two  possibilities  now 
exist:  either  there  are  two  independent  manifolds  corresponding  to  X,  or  there  is  only  one. 

If  two  independent  manifolds  exist,  then  they  span  the  plane  and  so  every  vector  forms 
a  manifold  with  this  same  eigenvalue  X.  If  X±0,  all  trajectories  are  straight  lines  through  the 
corresponding  equilibrium  position  which  is  now  termed  a  "star",  stable  or  unstable  if  A<0 
or  if  X>0,  respectively: 


Real,  Equal  Eigenvalues: 
A,=A2=A>0  "Star" 


If  the  eigenspace  corresponding  to  the  eigenvalue  X  is  one  dimensional,  i.e.,  there  is 
only  one  manifold  corresponding  to  X,  the  corresponding  fixed  point  is  a  "degenerate  node": 
as  t-±°°  where  all  trajectories  become  parallel  to  the  one  available  eigendirection.  Essentially, 
this  behavior  results  when  two  distinct  manifolds  are  scissored  together  with  some  of  the 
trajectories  becoming  trapped  in  the  collapsing  region  between  the  two  eigendirections,  while 
the  surviving  trajectories  are  pulled  around  to  form  the  degenerate  node  with  stability 
determined  according  to  X<0  or  A>0.  The  degenerate  node  is  a  borderline  case  between  a 
spiral  and  a  node;  the  trajectories  are  trying  to  wind  around  in  a  spiral  (complex  eigenvalues), 
but  insufficient  energy  is  present  to  complete  the  winding: 
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Real,  Equal  Eigenvalues  with  One 
Manifold:  X12  =A 

5.         Miscellaneous 

a.  Small  Non-linear  Terms  Change  Linearization  Stability-Sometimes 

In  the  presence  of  non-linear  terms,  the  character  of  marginally  stable  fixed 
points  can  change:  centers,  stars  and  degenerate  nodes  all  satisfy  the  criterion  Re(A)=0.  For 
example,  if  the  linearization  predicts  a  center,  the  influence  of  non-linear  effects  causes  the 
phase  portrait  to  show  a  stable  or  unstable  spiral,  a  stable  or  unstable  node;  the  presence  of 
non-linear  terms  violate  Re(A)=0  by  either  adding  or  removing  energy  from  the  system, 
driving  it  away  from  the  conserved  state  (center).  Similarly,  stars  and  degenerate  nodes  can 
also  change  their  character,  but  unlike  centers,  non-linear  effects  change  only  the  character 
of  the  star  or  degenerate  node,  not  their  stability:  stable  star  to  stable  spiral  or  stable  node, 
for  example. 

The  phase  portrait  is  structurally  stable  if  its  topology  is  not  changed  by  arbitrary  small 
perturbations  to  the  vector  field:  the  phase  portrait  of  a  saddle  point  is  structurally  stable,  but 
that  of  a  center  is  not  since  an  arbitrary  small  amount  of  damping  converts  the  center  to  a 
spiral.  Essentially,  if  the  phase  portrait  changes  its  topological  structure  as  a  parameter  is 
varied,  a  bifurcation  is  in  progress. 

b.  Example:  A  Glider 

Consider  a  glider  flying  at  speed  v  at  an  angle  6  to  the  horizontal.  Its  motion 
is  governed  approximately  by  the  dimensionless  equations 
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dv 
dt 


=  -sin(9)-Dv: 


vdQ.  =  -cos(e)  +  v2 
dt 


where  the  trigonometric  terms  represent  the  effects  of  gravity  and  the  v2  terms  represent  the 
effects  of  drag  and  lift.  Here,  it  is  desired  to  find  the  effects  on  the  glider  performance  as  the 
drag  parameter,  D,  is  varied. 

To  begin  the  analysis,  first  a  local  analysis  is  performed  about  any  stationary  points 
(fixed  points)  that  occur  when  the  derivative  terms  are  set  to  zero.  These  fixed  points  are 
then  used  to  evaluate  the  linearization  (Jacobian)  so  that  the  system  characteristics  about  the 
stagnation  can  be  determined  in  terms  of  eigenvalues  and  eigenvectors.  The  results  are 
summarized: 


Fixed  Point 

Comments 

Jacobian  Eigenvalues 

(roU) 

n=0,2,4,6,... 

K2-           2 

(mr,-l) 

n=0,2,4,6,... 

.         2Z)±^/4Z)2-8 

K2-            2 

From  the  preceding,  it  is  seen  that  if  the  drag  D>2  in  non-dimensional  terms,  the  eigenvalues 
of  both  fixed  points  will  always  be  real  and  distinct.  The  associated  dynamical  behavior, 
locally,  can  therefore  be  characterized  as  a  stable  node  for  the  fixed  point  (n7r,l)  and  as  an 
unstable  node  for  (n7i,-l).  If  D>0,  then  the  trace  of  the  Jacobian  (sum  of  respective 
eigenvalues)  will  always  be  less  than  zero  for  the  fixed  point  (nrc,l),  while  the  instability  for 
(n7r,-l)  is  shown  by  the  exponential  growth  of  the  oscillation  component  of  the  associated 


92 


system  eigenvalues  (trace>0). 

If  2>D>0,  the  eigenvalues  are  complex  and  hence  the  local  dynamics  of  the  system 
is  represented  by  stable  spirals  for  (nir,l)  and  an  unstable  spiral  for  (n7r,-l).  Similarly,  the 
stability/instability  is  evident  through  the  exponential  decay/growth  of  the  oscillations 
characterized  by  the  real  component  of  the  respective  eigenvalues.  It  should  be  noted  that 
when  D=2,  the  eigenvalues  are  real  and  equal,  however.  Consequently,  the  system  behavior 
should  be  that  of  a  star  or  a  degenerate  node:  a  spiral  (unstable  or  stable,  depending  on  the 
eigenvalues),  is  found  instead.  Here  is  evidence  that  nonlinear  terms  in  the  system  can  change 
the  local  behavior  from  that  which  is  predicted  by  the  linearization  analysis. 

If  D=0,  the  case  of  no  drag,  then  the  eigenvalues  for  the  fixed  point  (n7i,l)  with  n=0 
will  be  imaginary,  reflecting  pure  oscillation  of  the  system  and  consequently  periodic  solutions 
with  period  T=27r/a>.  The  resulting  equilibria  are  characterized  as  centers  and  describe  the 
lowest  energy  or  power  state  of  the  system,  the  glider  moving  in  a  potential  field,  gravity, 
without  any  other  influences.  That  is,  the  glider  moves  along  lines  of  constant  energy,  along 
a  curved  flight  path  with  a  radius  of  curvature  equal  to  the  absolute  altitude,  experiencing  the 
standard  "fictitious"  forces  due  to  relative  motion  (centripetal  acceleration,  for  example). 

In  general,  for  D>0,  the  case  of  increasing  drag,  the  centers  at  (n:r,l)  transition  to 
stable  spirals/nodes,  whereas  the  centers  at  (mr,-l)  transition  to  unstable  spirals/nodes.  From 
an  eigenvalue  consideration,  both  dampening  and  oscillations  are  now  evident;  the  eigenvalues 
are  complex  with  exponential  decay  for  the  stable  spiral/node  and  exponential  growth  for  the 
unstable  spiral/node.  In  essence,  drag  tends  to  offset  the  lift  generated,  as  expected.  This 
example  shows  that  if  the  phase  portrait  changes  its  topological  structure  as  a  parameter  is 
varied,  a  bifurcation  is  in  progress. 
D.         BIFURCATIONS 

As  indicated,  the  stability  of  a  fixed  point  can  change;  such  qualitative  changes  in 
dynamics  are  referred  to  as  bifurcations.  Bifurcations  provide  models  of  transitions  and 
instabilities  as  some  control  parameter  is  varied,  the  onset  of  coherent  radiation  in  a  laser,  for 
example.  Of  the  many  types  of  bifurcations  possible,  only  Saddle-node  and  Transcritical,  are 
outlined  here. 
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The  saddle  and  transcritical  bifurcations  occur  when  one  of  the  eigenvalues  equals 
zero,  involving  the  collision  of  two  or  more  fixed  points.  These  types  of  bifurcations  have 
analogs  in  one,  two  and  higher  dimensions  since  the  important  behavior  for  these  dynamics 
are  confined  to  a  one-dimensional  subspace  along  which  the  bifurcation  occurs,  when 
Re(A)=0,  while  in  the  extra  dimensions  the  flow  is  either  simple  attraction  or  repulsion  from 
that  subspace. 

1.  Saddle-Node  Bifurcations 

The  saddle-node  bifurcation  is  the  basic  mechanism  by  which  fixed  points  are  created 
and  destroyed.  For  instance,  as  a  parameter  is  varied,  two  fixed  points  move  toward  each 
other,  collide  and  mutually  annihilate.  The  prototypical  example  of  a  saddle-node  bifurcation 


.  7 

■y*  —  v^y  * 

is  given  by  the  system         .  _  _  where  r  is  a  small  parameter  determined  by  the 


problem  at  hand.  Here,  the  motion  is  decoupled  with  the  y-direction  assumed  arbitrarily  as 
exponentially  decaying.     The  fixed  points  for  this  system  are       (x*9y*)=  (±\/r,0)     , 


stagnation  points  which  exist  in  the  x,y  plane  when  r>0,  coalesce  when  r=0  and  do  not  exist 
(in  x,y  space)  when  r<0.  Upon  further  analysis,  it  is  seen  that  the  bifurcation  is  fundamentally 
one-dimensional,  with  the  fixed  points  sliding  toward  each  other  along  the  unstable  manifold 

of  the  saddle  point  at      (x  *y  *)  =  (--y/^O)     where  r>0. 

For  two-dimensional  system,  and  for  higher-dimensional  systems,  the  flow  is  limited 
not  only  by  fixed  points  (stagnation),  buts  also  by  closed  orbits  and  the  unions  of  fixed  points 
and  the  trajectories  connecting  them.  The  latter  are  referred  to  as  heteroclinic  orbits  when 


94 


they  connect  distinct  points  and  homoclinic  orbits  when  they  connect  a  point  to  itself.  For 
example,  if  the  physical  system  is  such  that  the  unstable  manifold  of  a  saddle  node  intersects 
with  its  stable  manifold,  the  resulting  orbit  is  called  "homoclinic:"  the  connection  of  the 
unstable  manifold  (X>0)  and  the  stable  manifold  (A<0)  of  the  saddle  point  is  common  in 
conservative  systems  but  can  also  occur  in  non-conservative  systems.  When  the  system 
dynamics  are  such  that  additional  damping  is  introduced,  the  tendency  is  for  the  homoclinic 
orbit  to  break,  causing  the  unstable  manifold  to  have  components  which  approach  the  fixed 
point  at  the  origin  of  the  phase  portrait  as  t-°°.  This  fixed  point  is  then  a  sink,  with 
eigenvalues  -a±iA  where  a  is  the  damping.  [Ref.  55:  p.  45] 

2.  Transcritical  Bifurcation 

There  are  certain  scientific  situations  where  a  fixed  point  (stagnation)  must  exist  for 
all  values  of  a  parameter  and  can  never  be  destroyed.  For  example,  in  the  logistic  equation 
there  is  a  fixed  point  at  zero  population,  regardless  of  the  value  of  the  growth  rate;  the  laser 
rate  equation  (classical)  is  another  example  for  which  a  fixed  point  always  exists.  However, 
such  a  fixed  point  may  change  its  stability  as  the  parameter  is  varied.  The  transcritical 
bifurcation  is  the  standard  mechanism  for  such  changes  in  stability,  with  normal  form  given 

Y  —  f  Y~^~  Y 

by:  where  r  is  control  parameter  determined  by  the  problem.    The 

y=-y 


difference  between  the  saddle-node  and  transcritical  bifurcations  is  that  in  the  latter,  the  fixed 
points  do  not  "disappear,"  their  stability  changes  instead.  Consider,  for  example,  the  laser 
rate  equation  [Ref.  54:  p.  54]: 

—  =gain  -  loss-GnW-kn 
dt 

- Gn{Wo-an) - kn  =  GWn  -  aGn2-kn 
=  (GWo-k)n-(aG)n2 

where  n(t)  represents  the  number  of  photons,  W  characterizes  the  rate  with  which  a  single 
excited  atom  spontaneously  generates  a  photon  per  second  offset  by  the  number  of  photons 
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per  second  that  drop  back  to  the  ground  state,  a;  W0  is  the  pump  strength.  G  is  the  gain 
coefficient,  and  -k  shows  the  rate  at  which  photons  are  lost  in  the  laser  by  scattering  or 
impurities.      For  the  laser  rate   equation   shown,   two   equilibria   exist:    n,*=0   and 


n2  = 


GWQ-k 
aG 


.  But  what  does  this  mean?   Shown  below  are  the  phase  portraits  with 


the  fixed  points  indicated: 


# 


-n 


/ 


/~ 


-n 


W<k?G 


w=i:'C 


w,  >k/a 


When  W0<  k/G,  GWo-k<0,  ni*=0  is  stable  and  for  low  levels  of  incident  energy  (weak 
pumping)  the  solution  %*=<)  must  be  valid  since  n(t),  the  number  of  photons  generated,  must 
always  be  positive.  Consequently,  there  can  be  no  laser  action  -  the  laser  acts  like  a  lamp. 
When  the  laser  is  pumped  at  a  greater  intensity  such  that  now  (GWo-k)>0  or  W0>k/G  holds, 
the  solution  with  n>0  becomes  possible  and  laser  action  is  attained.  The  system  undergoes 
a  transcritical  bifurcation  when  W0=k/G,  the  laser  threshold  for  this  model.  When  W0>k/G 
the  lamp  becomes  a  laser  and  the  fixed  point  n^O  loses  stability,  driving  the  system  toward 
n2*.  The  pumping  intensity  elevates  the  inversion  population:  a  critical  threshold  exists  when 
W0=k/G  below  which  laser  action  cannot  occur  (insufficient  pumping  energy  and 
consequently  insufficient  inversion)  and  the  laser  acts  like  a  lamp.  Above  the  threshold,  laser 
action  occurs  since  sufficient  energy  is  now  available  to  raise  the  atoms  to  higher  energy 
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levels  so  that  the  number  of  photons  generated  are  enough  to  cause  lasing.  At  the  threshold, 
at  the  transcritical  bifurcation  point,  the  stability  of  the  system  origin  becomes  unstable, 
driving  the  dynamics  toward  greater  population  inversion  and  consequently  toward  the  second 
fixed  point.  In  fact,  there  are  certain  scientific  situations  where  a  fixed  point  (stagnation)  must 
exist  for  all  values  of  a  parameter  and  can  never  be  destroyed  although  it  changes  its  stability 
character;  the  laser  is  such  an  example,  for  it  is  the  stagnation  at  the  origin  that  drives  the 
system  toward  the  lasing  threshold  at  the  bifurcation  point.  [Ref.  63] 
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APPENDIX  D.  CHARACTERISTICS  OF  SHEATH  EQUATIONS 


The  Jacobian  to  equations  (7-9), 
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has  the  following  characteristic  polynomial: 


characteristic   polynomial  =  aX5  +  pA.4  +  yX3   +  £A2  +  r\X  +  e  =  0 


where 
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To  find  the  solutions  to  the  characteristic  polynomial,  at  least  one  of  the  roots  has  to 
be  known  apriori.   To  this  end,  the  Jacobian  is  evaluated  for  the  fixed  point  found  when 

ne'  =n'  =n*  =  0     (see  Table  5  in  the  main  body  of  this  work)  but  with     E    retained  as 

a  non-zero  variable  (see  Chapter  V).  The  resulting  Jacobian  eigenvalues  for  the  fixed  point 
for  this  condition  become: 
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A,=[a]£         i2>3  =  0 


x.  ,  =  -l£H±N 


4,5 


2D. 


To   find    the   general   eigenvalues   of  the   Jacobian   then,   the   characteristics 
Aj^tfJ-E"    and    X2  =  0     are  assumed  as  solutions.    By  factoring  the  characteristic 

polynomial  with  these  results  and  neglecting  any  remainder  terms  gives  the  approximate  and 
complicated  solutions  to  the  characteristic  polynomial: 


E     A2  =  0 


+  2 


with 


101 


i   ( p    )     i        1 1  p    V 

6    %       J        2    2       27V  a       / 

i       i/p    v 

v=— Yi  _  —  — +K 

3    '       9l  a       / 


where 


K  = 
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If  numeric  data  nitrogen  is  applied  to  the  stagnation  region  near  the  plasma  boundary  (Tables 
4,  5)  and  the  results  substituted  into  these  roots,  an  order  of  magnitude  comparison  can  be 
performed  with  the  actual  results  obtained  (the  exact  evaluation  of  Jacobian  at  the  fixed  point 
in  question): 


K 

A2 

K 

\ 

K 

Jacobian  Evaluated  at 
Plasma  Boundary 
Stagnation 

1.376704 

0 

0.630455 

-1.37437 

-0.63279 

Approximate  Solutions 
to  Jacobian  at  Plasma 
Boundary  Stagnation 

1.41467 

0 

0.54 

-1.41262 

-0.54205 

Table  1 1 :  Sheath  Equations  -  Jacobian  Characteristics  Evaluated  at  Plasma  Boundary 

Equilibrium 

It  seems  that  Xl  dominates  the  exponential  growth  seen  during  any  numeric  simulation  -  the 
process  of  ion  repulsion  from  the  positive  plate  forces  the  instability  of  the  sheath. 
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